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ABSTRACT 

Chaos  metrics  are  examined  as  a  tool  to  analyze 
atmospheric  three-dimensional  dispersion  models  at  the 
individual  particle  rather  than  the  aggregate  level.  These 
include  the  self -af fine  fractal  dimension,  DA/  Shannon 
entropy,  S,  and  Lyapunov  exponent,  X.  Intercomparison  of  these 
metrics  is  first  performed  with  the  one-dimensional  logistics 
difference  and  the  two-dimensional  Henon  systems  of  equations. 
The  fractal  dimension  and  Shannon  entropy  are  then  measured  as 
a  function  of  the  inverse  Monin-Obukhov  length  (1/L)  for  two 
three-dimensional  Lagrangian  particle  dispersion  models,  the 
McNider  particle  dispersion  model  and  the  NPS  particle 
dispersion  model  now  under  development.  The  fractal  dimension 
and  Shannon  entropy  uncover  weaknesses  in  both  models  which 
are  not  obvious  with  standard  geophysical  measures.  They  also 
reveal  similarities  and  differences  between  the  atmospheric 
models  and  simple  chaos  systems.  Combined,  these  chaos 
measures  may  lend  detailed  insight  into  the  behavior  of 
Lagrangian  Monte  Carlo  dispersion  models  in  general. 
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I .  INTRODUCTION 

A.   OVERVIEW 

Unlike  the  relatively  simple  linear  relationships  studied 
in  basic  physics,  natural  phenomena  often  exhibit  complex 
nonlinear  behavior.  One  familiar  consequence  is  the  inability 
of  meteorologists  to  make  accurate  long-term  weather  forecasts 
(Devaney,  1990) .  Like  turbulence  in  other  systems,  atmospheric 
weather  systems  seem  to  follow  patterns  which  are  chaotic  and 
unpredictable  in  the  long  term.  The  behavior  of  a  particle 
released  in  such  a  system  is  deterministic  rather  than 
stochastic.  However,  small  differences  in  initial  conditions 
result  in  particle  paths  which  diverge  exponentially  with 
time.  Thus,  since  the  resolution  and  accuracy  of  measurements 
of  the  initial  conditions  are  always  limited,  forecasts  will 
diverge  from  reality  exponentially  with  time  such  that  long 
term  predictability  becomes  impossible  (Gleich,  1987) . 

In  fluids,  predictability  relates  to  particle  diffusion. 
Although  the  relationship  is  not  simple,  less  predictable 
systems  tend  to  more  diffusive. 

A  good  Lagrangian  particle  model  which  can  accurately 
simulate  diffusion  over  a  wide  range  of  atmospheric  conditions 
is  of  value  in  predicting  toxic  dispersion  from  various 
sources,  such  as  biochemical  warfare  agents,  Navy  LNG  tanker 


leaks,  nuclear  reactor  and  weapons  plumes,  industrial  chemical 
plants,  large  rocket  and  missile  launch  emissions,  aborts,  and 
static  test  firings,  and  a  variety  of  planned  and  emergency 
releases  from  tanks  and  hoses  involved  in  storage,  transfer, 
and  transport  of  dispersive  volatile  toxins.  Presently, 
predictions  using  existing  particle  models  may  vary  greatly 
from  reality  when  operating  in  certain  atmospheric  conditions. 
(Venkatram  and  Wyngaard,  1988) . 

Heretofore,  atmospheric  particle  dispersion  models  have 
generally  been  compared  against  turbulence  data  in  terms  of 
spectra  and  standard  geophysical  turbulence  measures,  such  as 
the  turbulent  kinetic  energy  (tke) ,  its  component  velocity 
variances  (ctu2,  av2,  and  aw2)  ,  and  the  Brunt-Vaisala  frequency 
(BVF) ,  an  atmospheric  stability  measure.  The  aggregate 
diffusion  of  many  Lagrangian  particles  has  also  been  compared 
with  the  results  of  Gaussian  plume  statistical  theories  as 
well  as  measurements  of  dispersing  clouds  (Venkatram  and 
Wyngaard,  1988) .  However,  complex  but  wholly  predictable 
periodic  behavior  such  as  ocean  tides  may  appear  quite  similar 
to  truly  chaotic  or  random  behavior  when  using  such 
techniques.  Another  potential  set  of  tools  in  analyzing 
particle  dispersion  model  performance  are  the  recently 
developed  chaos  metrics  such  as  fractal  dimension,  DA,  Shannon 
entropy,  S,  and  the  Lyapunov  exponent,  X  (Baker  and  Gollub, 
1990)  . 


B.  OBJECTIVE 

This  study  examines  the  self-affine  fractal  dimension,  DA, 
and  compares  it  to  the  Shannon  entropy,  S,  and  the  Lyapunov 
exponent,  X.  These  metrics  are  first  applied  to  the  l-D 
chaotic  system  known  as  the  logistics  difference  relation  and 
the  chaotic  2-D  Henon  system.  These  simple  well  characterized 
systems  are  studied  to  determine  similarities  and  differences 
between  the  above  three  chaos  metrics.  Two  3-D  Monte  Carlo 
Lagrangian  scattering  routines  designed  to  simulate 
atmospheric  dispersion  are  then  studied  as  representative  of 
more  complex  real  world  systems. 

C.  WHY  USE  CHAOS  METRICS 

There  are  at  least  two  good  reasons  to  try  chaos  metrics 
on  dispersion.  One,  chaos  metrics  may  provide  a  means  to 
discriminate  between  periodic  wave  behavior,  chaos,  and 
differing  degrees  of  turbulence.  Second,  chaos  metrics  may 
provide  some  additional  insight  into  the  behavior  of  particle 
dispersion  in  3-D  scattering  routines. 

Kamada  and  DeCaria  (1992)  have  shown  that  though  nocturnal 
periodic  gravity  waves  have  quite  different  dispersion 
characteristics  from  turbulence,  the  two  cannot  be 
distinguished  readily  by  using  standard  atmospheric  turbulence 
measures  such  as  the  Brunt-Vaisala  frequency  (BVF) ,  vertical 
velocity  variance  (aj)  ,      turbulent  kinetic  energy   (tke) , 


buoyancy  length  scale  (lb)  ,  or  spectral  analysis  using  FFTs. 
The  self-affine  fractal  dimension,  DA,  was  shown  to  be  the 
only  facile  wave/turbulence  discriminant  tested.  Since  the 
Shannon  entropy  and  Lyapunov  exponent  are  two  other  standard 
chaos  measures,  they  might  also  be  tested  as  possibly  useful 
dispersion  measures. 

In  the  3-D  Monte  Carlo  scattering  routines  which  model 
atmospheric  dispersion,  such  as  the  McNider  Particle 
Dispersion  Model,  chaos  metrics  might  quickly  determine 
certain  model  characteristics  for  a  given  range  of  parameters. 
For  example,  in  examining  an  expanding  cloud  of  particles,  the 
Lyapunov  exponent  would  measure  the  divergence  rate  of  the 
particles.  The  Shannon  entropy  would  measure  the  evenness  of 
distribution  of  particles  within  the  expanding  cloud.  The 
self-affine  fractal  dimension  could  measure  the  total  apparent 
distance  a  particle  travels  in  a  time,  T,  as  a  function  of  the 
time  resolution,  e,  used  within  T.  This  might  give  an 
indication  of  the  range  of  scales  over  which  mixing  occurs. 


II.   THEORY 

A.  OVERVIEW 

There  are  several  methods  of  measuring  chaos  and 
turbulence.  Of  interest  for  the  present  purposes  are  the  self- 
affine  fractal  dimension,  DA,  the  Shannon  entropy,  S,  and  the 
Lyapunov  exponent,  X.  Some  standard  geophysical  measures  will 
also  be  used  for  comparison  purposes  in  the  3-D  case. 
Descriptions  of  each  are  given  below  and  the  logistics 
difference,  Henon,  and  atmospheric  particle  diffusion  systems 
to  which  they  are  applied  are  then  described. 

B.  FRACTAL  DIMENSION  (DA) 

The  self-similar  fractal  dimension,  DB,  can  be  described 
readily  with  the  Cantor  set  (Devaney,  1990)  .  To  form  this  set, 
a  straight  line  is  drawn  with  the  middle  third  removed,  as  in 
the  second  line  of  Figure  1.  The  two  resulting  straight  lines, 
which  are  one-third  the  original  line  length,  are  then 
similarly  subdivided.  The  four  new  lines,  all  l/9th  the 
original  line  length,  are  further  subdivided,  ad  infinitum. 
From  top  to  bottom  in  Figure  1,  the  resulting  sets  of  lines 
are  self -similar ;  that  is,  the  manner  in  which  the  geometry  is 
altered  is  repeated  at  all  successive  levels  of  resolution. 

The  Koch  snowflake  is  another  geometrically  self-similar 
figure  (Devaney,  1990) .  It  is  constructed  by  removing  the 


'igure  1,  Construction  of  the  Cantor  Set 


middle  third  of  each  side  of  an  equilateral  triangle,  and 
replacing  it  with  two  pieces  of  equal  length,  creating  a  six- 
pointed  star.  This  star  has  12  sides,  all  of  length  1/3.  The 
middle  third  of  each  of  the  12  sides  is  again  removed,  and 
replaced  with  two  lines  of  length  1/9.  Again,  the  process  is 
repeated  ad  infinitum. 

Like  the  Cantor  set,  the  Koch  snowflake  is  also  self- 
similar  [Figure  2].  The  jagged  sides  appear  geometrically 
similar  at  increasing  levels  of  resolution. 

The  self-similar  fractal  dimension  is  defined  by  (Devaney, 
1990)  as 

D         In  (Total  length  of  pieces)  fli-i} 

B  In  (resolution) 

In  the  Cantor  set  example,  the  total  length  of  segments  of 
unit  length  at  level  n  is  2n,  and  the  resolution  level  is  3n, 
so  that 

=  _ln_lf  =  n  }n  2    *   0.6309  .  (H-2) 

B       In  3n   n   In  3 

Similarly  for  the  Koch  snowflake,  there  are  4  pieces  for 

each  level  of  n  with  a  magnification  of  3 ,  so  that 

=  In  An   m  lm262   .  (II-3) 

B        In  3n 

With  geometric  figures,   DB  is  unitless;  the  measure 

involves  a  length  divided  by  a  length.  However  this  definition 

of  fractal  dimension  cannot  apply  directly  to  a  time  series 


Figure  ,.  construction  of  the  Koch  Snowflafce.  Pro,  Devaney 

(1990)  . 


trace,  since  it  would  involve  the  square  root  of  the  amplitude 
squared  plus  the  time  squared,  i.e., 


In 


D„   = 


£  JAA2   +  At: 


(II-4) 


3  In  fit 

This  definition  must  be  adjusted  when  applied  to  time  series 
data  to  ensure  that  the  end  result  is  independent  of  the 
arbitrary  unit  scaling  between  amplitude  and  time. 

There  are  other  ways  to  characterize  the  fractal  dimension 
of  a  system.  For  a  more  suitable  times  series  measure,  McHardy 
and  Czerny  (1987)  redefined  the  length  metric  as 


L(e)   =  If  |F(t  +  e)  -  F(t)\  dt   .  (H-5) 

e  J 


G 

0 

where  F  is  the  amplitude  of  the  time  series  at  time,  t,  e  is 
the  time  increment,  and  T  is  the  time  window  over  which  L  is 
defined. 

This  definition  effectively  avoids  the  units  scaling 
problem.  Since  the  inverse  time  in  1/e  cancels  the  time  units 
in  the  integral,  L(e)  is  only  in  units  of  amplitude.  The 
fractal  dimension  is  then  defined  as 

n  _    din  Lie)  (II_6) 

.      dine 

McHardy  and  Czerny  applied  these  definitions  to  the  time 
variance  of  X-ray  luminosity  data  from  the  Seyfert  galaxy 
NG5506.  Collins  and  Kastogi  (1989)  later  applied  McHardy  and 


Czerny's  definitions  in  analyzing  gravity  wave  spectra  from 
the  atmospheric  mid-troposphere.  Recently,  Kamada  and  DeCaria 
(1992)  applied  DA  as  a  tool  for  discriminating  waves  from 
turbulence  in  nocturnal  atmospheric  boundary  layer  data. 

To  actually  compute  this  function,   the  integral  is 
approximated  with  a  numerical  summation,  so  that 

T 

L(e)   =  ££  £|F(t-e)  -  F(t)\    .         (Ii-V) 

G   o 

Since  at  each  resolution  St=e ,  the  leading  term  on  the  right- 
hand  side  is  always  unity,  so  that  the  length  is  now 


Lie)    =  £|F(t-e)  -  F(t)\     .  (II-8) 


Thus  L(e)  is  the  total  amplitude  change  over  a  time  series  of 
length,  T,  for  a  given  time  resolution,  e.  In  this  form  it  is 
quite  clear  that  time  is  removed  from  the  length 
determination,  so  that  L(e)  does  not  depend  on  some  arbitrary 
scaling  between  amplitude  and  time. 

C.   SHANNON  ENTROPY  (S) 
1.   Definition 

When  a  particle  is  first  released,  its  initial 
location  is  known  and  completely  specified,  so  the  information 
entropy  (defined  later)  for  its  location  is  zero.  Later, 
according  to  given  equations  of  motion,  its  position  diverges 
from  the  initial  point.  If  the  range  of  its  possible  positions 
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is  partitioned  into  equal  segments,  then  for  a  given  time 
interval,  T,  its  motion  can  be  recorded  in  terms  of  occupancy 
time  for  each  segment.  Then  the  particle  dispersion  rate  might 
be  measured  by  the  seeming  degree  of  randomness  of  occupancy 
time  or  evenness  of  the  state  probability  distribution  over 
time  period,  T.  Shannon  or  information  entropy  is  defined  by 
the  state  probability  distribution,  so  entropy  may  be  regarded 
as  a  measure  of  dispersion  rate  for  a  time  series  of  particle 
states.  This  can  apply  to  the  actual  particle  position,  its 
velocity,  or  its  phase  velocity.  Shannon  entropy  is  defined 
simply  as  (Baker  and  Gollub,  1990) 


N 

I 

i=l 


S  =  -  £  Pi  lnPi  ,  (H-9) 


where 


S  =  system  entropy, 
N  s  number  of  permitted  states,   and 
Pj  =  probability  of  state  i,   such  that  J^  pi  =  1.0 


N 


1=1 


Then,  for  N  permitted  states  (or  position  intervals) , 
the  maximum  possible  entropy  corresponds  to  equal  occupancy 
time  in  each  of  the  N  states,  so  p;  is  1/N  for  all  i.  That  is, 

S_  =  -  T    4  In  4  •  (II-10a) 


'max 


=  -  E  -  ln  -   ■ 

pi  N  N 


Expanding  the  summation  results  in 
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0   _  1 ,  1   l1nl   l1nl 

■SL=.„  =  — In —  +  — In —  +  — In  — 
max    N        N        N        N        N        N 


and  collapsing  the  common  multiples  gives 


(II-10b) 


'max 


=  -iWA  ml) 

\N         Nf 


(II-lOc) 


Smax  "  In  N 


which  corresponds  to  total  randomness,  or  a  completely  even 
particle  distribution  across  all  allowable  states.  Also, 


min 


=  o  , 


(11-11) 


which  corresponds  to  all  particles  being  in  one  state. 

For  a  randomly  moving  particle  in  2-D  space,  the 
computation  is  similar.  The  2-D  space  may  be  divided  into, 
say,  a  100  X  100  grid.  If  the  particle  is  moving  completely 
randomly,  at  time,  t,  it  may  be  in  any  one  of  the  grid  boxes 
with  equal  probability.  The  entropy  for  this  system  is  then 


S  =  - 


1002 


In 


100' 


5  =  -  1002 


100' 


In 


1002 

\ 


In 


ioo2; 


1002J 


which  is  the  same  as 


(II-12a) 


S  =   In  1002  =  St 


(II-12b) 
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2 .   Shannon  Entropy  for  an  N-D  System 

For  a  1-D  system,  the  domain  is  simply  partitioned 
into  n  intervals,  and  the  probability  computed  for  each 
interval.  For  a  simple  system  such  as  recursion  of  the 
logistics  difference  equation, 

xn*X   =  |*Xa(l-Xa)  .  (H-13) 

the  probability  for  the  ith  interval,  p;,  is  simply  the  number 
of  times  the  interval  has  been  occupied  after  n  number  of 
recursion  steps,  divided  by  the  total  number  of  steps. 

Note  that  the  number  of  partitions  should  be 
appropriate  for  the  total  number  of  steps.  Having  too  few 
intervals  is  equivalent  to  too  low  a  resolution  of  the  entropy 
and  may  result  in  a  misleadingly  low  value  of  S.  For  instance, 
with  only  one  interval,  pj  =  1,  and  S  =  0.  Ideally,  the  number 
of  intervals  for  maximum  resolution  of  the  entropy  of  the 
system  is  e^,  where  X  is  the  positive  Lyapunov  exponent  for 
the  system  and  N  is  the  number  of  steps  (Baker  and  Gollub, 
1990,  pp. 126-129) .  Since  X  is  typically  of  order  unity,  e^  can 
easily  be  a  computationally  intractable  number. 

Again,  for  3-D  systems,  assuming  N  equal  segments 
along  each  axis,  the  number  of  partitions  will  be  N3,  so 
computation  quickly  becomes  unwieldy  for  large  N.  In  practice, 
N  of  order  102  to  103  has  been  used  by  some  authors  (Baker  and 
Gollub,  1990,  pg.  88)  as  a  compromise  between  entropy 
resolution  and  computational  efficiency.  In  analyzing  real 
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data,  the  number  of  partitions  should  be  related  in  some 
direct  fashion  to  the  maximum  resolution,  e{,  of  the  measuring 
devices,  and  the  total  measurement  time,  T.  The  total  number 
of  partitions  should  probably  not  exceed  T/eif  or  otherwise 
even  a  completely  random  distribution  would  still  result  in 
some  intervals  unoccupied. 
3  .   Ln  vs  Log, 

Wolf  maintains  that  the  Shannon  entropy  should  be 
computed  with  log2  rather  than  the  natural  logarithm,  loge 
(Wolf,  1986,  pg.  276).  With  log2/  the  Shannon  entropy  relates 
directly  to  information  in  the  form  of  binary  bits,  since  one 
bit  of  information  can  be  specified  as  being  in  only  one  of 
two  states,  e.g.,  true  or  false,  on  or  off,  one  or  zero.  That 
is,  the  log2  basis  sets  the  entropy  equal  to  the  minimal 
length  of  binary  code  required  to  describe  the  state  of  the 
system.  If  all  particles  occupy  only  one  state,  turning  the 
bit  for  that  state  to  the  "on"  position  is  sufficient  to 
specify  the  state  of  the  system.  If  the  particles  are  evenly 
distributed  among  all.  states,  the  length  of  binary  code 
required  to  specify  the  system  is  equal  to  the  number  of 
states,  which  means  that  the  system  is  completely  random 
(Tribus  and  Mclrvine,  1971) .  However,  for  the  purposes  of  this 
study,  the  Shannon  entropy  can  also  be  written  in  the 
following  form: 
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N 


s  =  -k"£  Pilnpt    .  (H-14) 


2-1 


The  only  difference  in  measuring  entropy  by  the 
natural  logarithm  versus  base  two  is  the  value  imposed  on  the 
constant  K.  Most  computations  assume  that  K=l.  As  long  as  the 
method  of  computing  S  is  similar  when  making  comparisons  of 
computed  entropy,  there  need  be  no  confusion. 

D.   LYAPUNOV  EXPONENT  (X) 

For  an  expanding  turbulent  cloud  of  particles,  one  measure 
to  describe  the  system  is  the  particle  divergence  rate.  Due  to 
turbulence,  the  trajectories  of  two  particles  arbitrarily 
close  will  diverge  with  time.  The  divergence  rate  can  be 
characterized  by  Lyapunov  exponents.  The  Lyapunov  exponent,  X, 
is  also  a  measure  of  the  system's  sensitivity  to  initial 
conditions.  In  an  N-dimensional  system,  there  will  be  N 
Lyapunov  exponents;  however,  these  do  not  necessarily 
correspond  to  coordinate  axes.  So  in  a  Cartesian  coordinate 
system,  X,,  X2,  and  X3  are  locally  defined  Lyapunov  exponents 
which  generally  cannot  be  described  as  Xx,  Xy,  and  Xz. 
Direction  is  adjusted  for  each  point  along  the  trajectory.  In 
one-dimension,  the  Lyapunov  exponent  is  defined  by  (Wolf, 
1986,  p.  275)  as 
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JV-1 


X   =  lim  1  Y  In  If  (i)  |  ,  (11-15) 

which  can  be  described  over  a  map  domain  as  an  integral, 

N 

A.  =  f  pd   In  \f'(i)  |  di     .  (11-16) 

o 

Overall,  the  Lyapunov  exponent  can  be  viewed  as  a  measure  of 
the  average  local  stretching  rate  of  the  particle  trajectory, 
as  indicated  by  the  log  of  the  length  of  the  slope,  |f'(i) |, 
weighted  by  the  probability  of  encountering  that  slope. 
(Wolf,  1986) 

The  Lyapunov  exponent  is  related  to  the  entropy,  S,  as 
well  as  the  information  loss  rate.  For  instance,  if  a  measured 
position  is  presently  known  to  a  precision  of  16  bits,  and 
X  =  2  bits  per  second,  then  the  particle xs  future  trajectory 
cannot  be  predicted  to  any  degree  of  precision  beyond  8 
seconds  (16  orbits). 

It  should  be  noted  that  the  Lyapunov  exponent  defines  the 
average  rate  of  loss  of  predictive  power,  and  may  vary  locally 
along  the  orbit. 

The  Lyapunov  exponent  is  readily  computed  for  one 
dimension,  but  in  higher  dimensions  calculation  becomes  more 
complex.  In  three  dimensions,  the  directions  of  trajectory 
divergence  and  contraction  must  be  defined  in  terms  of  local 
tangents,  which  vary  from  point  to  point,  so  the  calculation 
of  the  exponents  must  constantly  adjust  for  each  change  in 
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direction.  Wolf  has  developed  algorithms  for  computing  higher 
dimensional  Lyapunov  exponents  which  involve  reorienting  the 
major  axis  of  an  ellipse  for  each  point,  then  renormalizing 
the  function  after  every  few  points  so  that  the  unit  ellipsoid 
does  not  overlap  an  attractor  (Wolf,  1986) . 

In  3-D,  three  Lyapunov  exponents  are  required  to  classify 
the  system.  A  negative  exponent  indicates  a  dissipative 
dimension.  If  all  three  exponents  are  negative,  the  system  is 
dissipative,  e.g.,  a  pendulum  settling  down  to  a  fixed  point. 
If  an  exponent  is  zero,  the  system  orbits  about  a  fixed  point. 
If  one  of  the  three  exponents  is  positive,  the  system  is 
chaotic;  the  orbital  trajectories  are  diverging. 

E.   GEOPHYSICAL  MEASURES  OF  TURBULENCE 

1 .  Overview 

Geophysicists  use  several  standard  methods  to  describe 
turbulence.  Discussed  here  are  those  methods  used  to  examine 
the  McNider  and  NPS  dispersion  models. 

2.  Turbulent  Kinetic  Energy  (tke) 

The  mean  turbulent  kinetic  energy,  or  TKE,  is  defined 
as 


tke   =  -|(oJ  +  a\  +  ol)      ,  (H-17) 


where 


au  =  standard  deviation,  u  component  of  particle  velocity, 
av  =   standard  deviation,  v  component  of  particle  velocity, 
aw  =  standard  deviation,  w  component  of  particle  velocity. 
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3.   Gradient  Richardson  Number  (Ri) 

The  gradient  Richardson  number  is  defined  by 


Ri  = 


ae 

s  g dz (11-18) 


jdu\2       (dv\2 
5zJ    I dz) 


where  the  overbars  signify  time  averages,  and 

g  =   earth^s  gravitational  acceleration, 

z  =   vertical  position  above  the  surface, 

9   =   potential  temperature,  i.e.  the  temperature 

normalized  for  adiabatic  expansion  with  pressure,  so 

6  =  t(JP-\^     ,  (11-19) 

where  Rd  is  the  gas  constant  for  dry  air,  and  Cp  is  the 
specific  heat  at  constant  pressure. 

Ri  is  used  in  place  of  the  Reynolds  number  as  a 
dynamic  indicator  of  turbulence  when  the  atmosphere  displays 
a  non-neutral  density  profile.  The  numerator  is  proportional 
to  the  buoyant  production  or  destruction  rate  of  tke, 
depending  on  negative  or  positive  sign,  respectively.  The 
denominator  is  proportional  to  the  shear  generation  rate  of 
tke,  which  is  nearly  always  positive.  Thus,  a  positive 
Richardson  number  indicates  a  stable  atmosphere  (increasing 
potential  temperature  with  height)  and  suppression  of 
turbulence.  Commonly,  when  the  buoyancy  destruction  rate  of 
tke  exceeds  1/4  the  shear  production  rate,  turbulence  is 
suppressed,  i.e.,  when  Ri  >  1/4  (Stull,  1988,  pg.  176). 
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4.   Brunt-Vaisala  Frequency  (BVF) 

The  Brunt-Vaisala  frequency,  BVF,  is  a  measure  of 
static  stability.  BVF  is  defined  by 


BVF  = 


g-H^Z     ,  (11-20) 


>  *v    9Z 


and  in  principal  qives  the  hiqhest  qravity  wave  frequency 
which  a  fluid  can  support.  It  is  undefined  for  neqative 
temperature  qradients,  i.e.,  an  unstable  atmosphere. 
(Sorbjan,  1989,  pq.  35). 

5.  Buoyancy  Length  (1B) 

The  buoyancy  length,  lb,  is  the  standard  deviation  of 
the  vertical  velocity  divided  by  the  Brunt-Vaisala  frequency, 

1B   =  -^   .  (H-21) 

B        BVF 

The  buoyancy  length  is  meant  to  be  a  measure  of  the  dominant 
eddy  scale.  (Stull,  1988,  pg.  310) . 

6.  Monin-Obukhov  Length  (L) 

The  Monin-Obukhov   length   is   given   by 

L  = =—      ,  (11-22) 

gkwWQ 

where 
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w'uQ  =  surface  temperature  flux  , 
k  =  von  Karman  cons  tan  t  ~  0 . 4 


,  V 

and    u,  = 


,      z         .  (11-23) 

u.  is  the  friction  velocity,  a  measure  of  surface  drag 
due  to  turbulent  friction.  The  Businger  function,  \pm,  is  a 
stability  correction  which  is  approximated  with  a  rational 
function  by  (Kamada,  1992b) 

-  5—  ,  —  >   0  (stable)  , 

Li  -Li 


*-  = 


0.15961583  -  5.4151107  — 


/   \  2      T 

1  -  3.59002232—  -  0.799168457' 
L 


<.   0  [unstable)   . 


(11-24) 
In  the  convective  boundary  layer,  L  is  proportional  to 

the  height  at  which  the  buoyant  generation  rate  of  tke  matches 

the  shear  generation  rate.  L  is  the  primary  stability  measure 

for  the  atmospheric  surface  layer  and  is  also  used  in 

combination  with  the  inversion  height,  z;,  to  characterize 

boundary  layer  stability.  L  >  0  indicates  a  stable  boundary 

layer  where  vertical  turbulence  is  suppressed  by  positive 

density  gradients  with  height,  z.  L  <  0  indicates  an  unstable 

surface  layer  where  upward  vertical  motion  is  encouraged  by 

negative  density  gradients.  A  stable  atmosphere  implies  that 
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a  vertically  displaced  air  parcel  will  tend  to  return  to  its 
original  height.  (Businger,  1973) . 

F.   1-D  AND  2-D  CHAOS  EQUATIONS.  LOGISTICS  &  HENON  FUNCTIONS 

1.  Overview 

The  initial  focus  of  this  study  is  to  compare  and  test 
some  simple  potential  dispersion  metrics.  The  fractal 
dimension,  DA/  the  Shannon  entropy,  S,  and  the  Lyapunov 
exponent,  X,  are  simple  expressions  which  can  monitor 
transitions  between  periodic  and  chaotic  behavior,  akin  to  the 
transition  between  laminar  and  turbulent  behavior  in  a  real 
fluid.  Two  such  expressions,  well  characterized  in  the 
literature,  are  the  logistics  difference  equation: 

xn+1  =  \ixn(l-xn)      ,  (11-25) 

and  the  Henon  system: 

xn+i   =  1  "  axl  +  Yn     > 

(11-26) 

Yn+l   =  bxn      • 
(Gould  and  Tobochnik,  1988,  pp. 152-178;  Baker  and  Gollub, 
1990) . 

2.  The  Logistics  Difference  Equation 

Though  simple  and  one  dimensional,  recursive  iteration 
of  the  logistics  difference  equation  results  in  chaotic 
behavior  for  certain  parameters.  For  0  <  x0  <  1  and  ju  <  1,  xn 
converges  to  0.  For  1  <  /x  <  3 ,  and  the  same  range  for  x0,  xn 
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converges  to  /u/ 4.  For  3  <  jii  <  3.6,  x,,  fluctuates  between  2n 
discrete  points;  while  beyond  ji  —  3.6,  the  solution  is  nearly 
random,  i.e.,  chaotic.  [Figure  3] 

3.   Graphical   Analysis   of   the   Logistics   Difference 

Equation 

The  logistics  difference  equation  is  parabolic,  so 
that  for  an  initial  x0  (given  as  0.04  with  ju=2  . 9  in  Figure  4), 
drawing  a  vertical  line  to  the  parabola  corresponds  to  xl. 
Then  a  horizontal  line  drawn  to  the  inscribed  straight  line  of 
slope  unity  corresponds  to  recursing  x„+1  to  a  new  value  for 
xn.  Reiterating  this  procedure  marches  the  solutions  to  a 
final  value  of  0.655. 

This  final  value  can  be  determined  exactly  from  the 
original  function.  At  steady  state,  xn+1  =  xn.  Then  for  Figure 
4 ,  where  \i   =   2.9, 

xn+1  =  2.9xfl(l-xfl)  ,  (11-25) 

and  solving  for  x„: 

Xn  =  0  ,       Xn  =  0.65517      .  (11-26) 

Note  since  the  equation  is  quadratic,  that  there  are 
two  solutions.  The  x„=0  solution  is  not  stable  (as  defined 
below) .  The  general  rule  is:  if  the  magnitude  of  the  slope  of 
the  function  in  the  region  of  the  solution  is  greater  than 
45°,  the  fixed  point  is  an  unstable  solution;  if  the  slope  is 
less  than  45°,  the  fixed  point  is  a  stable  solution. 
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Figure  3.  Logistics  Difference  Equation:  1-D  Chaotic  Motion, 


^n  +  l 


=  3.72xn(l-xn)  ,  with  x0  =  0.01 
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Figure  4.  Logistics  Difference  Equation  Graphical  Analysis 
xn+1  =  2.9xn(l-xn)  . 
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Figure  5.  Logistics  Difference  Equation  Graphical  Analysis 
xn+1  =  3.9xn(l-xn) 


25 


When  the  function  is  chaotic,  there  are  still  only  two 
solutions  to  the  quadratic  equation.  However,  both  points  are 
unstable:  they  both  repel  rather  than  attract  the  value  of 
xn+1.  The  two  points  can  still  be  computed,  resulting  in 

xn   =  0,  xn   =  0.74   .  (11-27) 

Another  way  of  looking  at  this  is  by  taking  the  derivative 
of  xn+1  with  respect  to  xn: 

dXn+1    =  |i  (1  -  2xn)       .  (11-28) 


dxn 


Checking  the  chaotic  system  of  Figure  5,  when  x,j  =  0.74,  the 
slope  is 

dx 

=  3.9(1  -  2(0.74))  -  -1.87   ,       (11-29) 


dxn 


which  is  less  than  -1,  and  hence  unstable. 

This  concept  is  vital  when  visualizing  what  the  Lyapunov 
exponent  is  measuring.  Any  time  both  solutions  are  at  points 
where  the  slope  is  greater  than  45°  or  less  than  -45°,  the 
Lyapunov  exponent  is  greater  than  zero  (In  of  the  slope,  which 
is  greater  than  unity) ,  and  the  system  is  chaotic. 

4.   Bifurcation  Diagrams 

A  map  of  the  possible  values  of  xn  for  various  n  is 
called  a  bifurcation  diagram  [Figure  6] .  With  the  bifurcation 
diagram,  the  complete  behavior  of  the  function  can  be 
determined.  For  example,  for  0  <  /i  <  3 ,  xn  tends  to  only  one 
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Figure  6.  Logistics  Difference  Equation  Bifurcation  Diagram. 
xn+i  =  /xxn(l-xn)  . 


27 


value;  the  repetitive  recursion  of  xn+1  to  xn  points  to  one 
attracting  fixed  point,  corresponding  to  a  dissipative  system. 
At  ji  >  3 ,  the  bifurcation  diagram  splits  into  two  fixed 
points.  The  solution  of  xQ+1  bounces  back  and  forth  between 
these  two  points.  This  is  period  doubling  to  period  2.  Note 
for  still  larger  values  of  m  that  there  is  another  split  to 
period  4,  then  period  8,  followed  by  a  rapid  series  of 
bifurcations  culminating  in  chaos,  where  the  points  appear  to 
be  distributed  over  a  wide  range  of  xn. 

Another  curiosity  appears  in  the  bifurcation  diagram. 
There  are  numerous  "windows"  within  the  chaotic  region.  The 
most  obvious  is  around  p.  —  3.82,  where  the  period  4  and  8 
behavior  seems  to  appear  again.  Increased  resolution  would 
show  many  more  such  windows.  Also  note  that  period  widths  get 
progressively  narrower,  and  that  for  higher  periods  it  is 
difficult  to  distinguish  say  period  64  from  chaos. 
5.   The  Henon  Equations 

As  with  the  logistics  difference  equation,  the  Henon 
set  has  stable  solutions  for  a  range  of  parameters  a  and  b 
which  correspond  to  one  fixed  point.  Other  values  for  a  and  b 
result  in  two  fixed  points,  or  period  doubling  [Figure  7],  and 
period  4  and  higher.  For  still  higher  values  of  a  and  b,  the 
result  is  chaotic  behavior:  the  possible  positions  appear  to 
be  spread  throughout  a  discrete  range  of  x  and  y  [Figure  8]. 
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The  Henon  bifurcation  diagram  resembles  the  logistics 
difference  equation  diagram,  but  has  some  obvious  differences 
[Figure  9].  One,  it  does  not  show  symmetric  splitting,  but 
rather  is  skewed  downward  from  the  center line.  Two,  at  the 
value,  a  =  1.08  (with  b  =  0.3),  a  window  appears  in  the 
bifurcation  map  with  new  values  for  xn;  the  periodicity  in 
this  region  occurs  at  points  outside  the  previous  range  of  xn. 
These  new  points  seem  to  appear  out  of  nowhere,  and  are  not 
linked  to  previous  points  by  bifurcation. 

G.   3-D  ATMOSPHERIC  PARTICLE  DISPERSION  MODELS 
1.   Overview 

The  second  focus  of  this  study  is  to  examine  chaos 
metrics  as  performance  measures  for  3-D  Monte  Carlo  scattering 
routines.  Of  immediate  interest  is  particle  dispersion  within 
the  atmosphere.  The  McNider  Dispersion  Model  is  commonly  used 
to  estimate  the  transport  and  diffusion  of  atmospheric 
pollutants.  The  model  estimates  diffusion  based  on  atmospheric 
flow  predictions,  and  simulates  the  behavior  of  pollutant 
clouds  by  representing  them  as  the  ensemble  trajectories  of 
numerous  Lagrangian  point  particles  (Pielke,  1984) .  For  this 
study,  the  entire  range  of  possible  atmospheric  stabilities 
was  studied,  from  -0.2  <  1/L  0.1.  Therefore,  to  render  the 
computations  tractable,  the  flow  predictions  were  obtained 
from  boundary  similarity  theory  developed  by  Sorbjan  (1990) 
and  Kamada  (1992  a,b)  rather  than  from  a  mesoscale  windflow 
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Figure  9.     Henon  System:   Bifurcation  Diagram,   xn+1  =  l-axn2+yn, 
Yn+i   =   0-3xn 
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model,  since  steady  state  mean  flows  could  be  assumed.  The 
difference  is  not  large,  since  most  of  the  similarity 
algorithms  are  also  used  in  the  turbulence  closure  scheme  for 
the  windflow  model. 

Several  weaknesses  were  perceived  in  the  McNider 
formulation.  Therefore,  a  similar  model  was  developed  at  NPS, 
and  was  also  tested.  It  features  a  double  Gaussian  vertical 
velocity  skewness  scheme  and  damping  of  the  particle 
oscillations  within  the  boundary  layer  (Kamada  1992b) . 
2 .   The  McNider  Dispersion  Model 

For  this  study,  mean  flow  components  are  neglected  in 
order  to  focus  on  the  turbulent  fluctuations.  The  McNider 
Dispersion  Model  utilizes  the  form 

x(t  +  At)  =  x(t)    +  u(  t)  At, 

y(t  +  At)  =y(t)  +  v(t)At,         (11-30) 

z(t  +  At)  =  z(t)    +   v(t)At, 

where  x,  y,  and  z  define  the  particle  position  in  a  Cartesian 
coordinate  system,  and  the  u,  v,  and  w  terms  refer  to  the 
fluctuating  turbulent  velocity  components,  respectively.  The 
velocity  components  are  computed  by  (Pielke,  1984) 

u(t)  =  u(t  -  At)Ru(At)    +u/(t-At), 

v(t)    =  v(t  -  At)Rv(At)    +  v'(t   -  At)  ,      (H-31) 

w(t)    =  wit  -  At)RjAt)    +t/(t-At), 
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where  R  refers  to  the  Lagrangian  autocorrelation  factors  for 
each  velocity  component  and  is  a  function  of  At.  The  second 
terms  on  the  right  refer  to  random  fluctuating  components  of 
the  velocity,  and  are  related  to  the  flow  field's  turbulent 
kinetic  energy.  This  random  fluctuation  of  the  velocity  is 
constrained  by  the  fluid's  turbulent  kinetic  energy  with 


o„  = 


ou^l  -  i^(At)  , 
av]/l  -  Rl(bt)  , 
a'w  =  awSjl  -  Rl(At)  , 


ov  = 


(11-32) 


where  the  a  terms  refer  to  the  standard  deviations  of  the 
velocity  components.  These  parameters  are  not  directly 
computed  from  atmospheric  windflow  in  McNider's  model.  Instead 
an  empirical  formula  is  used  to  deduce  these  values: 


o     =  \ 

w 


where 


1.21 


(  Ric  -  Ri\ 


Ri, 


0.58 


i  an)2  +  (ixT 

\{dz)         { Bz) 


i*»  ■ 


i>°  ■ 


(11-33) 


1   = 


\0.35z, 
17  0  in, 


z  <  205m     , 
z  £  205fl]     . 


(11-34) 
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1^  is  the  local  exchange  coefficient,  described  later  in  the 
NPS  model  discussion.  To  account  for  advection  across  vertical 
gradients  of  ctw, 


°"  =  bt^iwU  ~ At)  +  °"k  ' 

and  the  horizontal  components  are  determined  by 


(11-35) 


ou  =  ov  =  \ 


u. 


12 


zi 


\2 

3 


2  L 


2.3u,  , 


*  S  0 
L 

4  >  0 


(11-36) 


The  maximum  wavelength  in  the  vertical  velocity  spectra  is 
given  by 


0  £  z<.   \L\     , 


**.  =  < 


0.55  +  0.38—  \ 
5.9  z  , 


1.8^ 


1  -  exp 


-4  2 


0.0003exp|  — 


L  <  z   £  0  . 1  z±     , 


0  . 1  zi  <  z  <  zi 


(II-37a) 


for  z/L  <  0,  and 


A.m  =  z  ;    Am  <;  2.9  1  ,   —  £  0   . 


(II-37b) 


In  the  above,  Z;  denotes  the  top  of  the  planetary  boundary 
layer. 

The  parameter,  Ri,  is  the  gradient  Richardson  number, 
defined  previously  by  eguation  11-19.  The  critical  Richardson 
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number,  Ric,  beyond  which  turbulence  is  suppressed,  is  given 
by 


Rir   =  0.115  Az0'175  , 


(11-38) 


where   Az   is   measured   in   centimeters.   The   Lagrangian 
autocorrelation  terms  are  determined  by 


Ru(At)    =  exp| 
Rv[At)    =  exp 

Rw(At)    =  exp 


-At" 
'  -At" 
f-At 


(11-39) 


TL  is  the  Lagrangian  integral  or  e-folding  time  scale  for 
velocity  autocorrelations,  and  is  determined  for  each 
component  from  the  turbulence  spectra  with 


(11-40) 


where  Xm  is  the  dominant  wavelength  for  each  component 
Furthermore,  the  average  velocity  is  defined  by 


V  =  Ju2  +  v2   +  w2    .  (11-41) 

13,    the  ratio  of  Lagrangian  to  Eulerian  time  scales,  is  given 
by 
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(11-42) 


where  /5  is  restricted  to  (3   <   10. 

The  horizontal  components  for  peak  wavelength  are 
given  by 


(11-43) 


The  constant,  A,  varies  as  a  function  of  stability  in  the 
unstable  boundary  layer,  and  is  computed  from 
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Values  of  u'(t-At)  and  v'(t-At)  are  determined  randomly  from 
a  normal  probability  distribution  with  zero  mean.  However,  the 
turbulence  distribution  of  vertical  velocity  is  skewed  [Figure 
10]  ,  and  the  normal  distribution  is  modified  by  changing  the 
method  of  computing  w: 
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(11-45) 


where 


a  =  -0.028  +0.6 \P\ 
r\   =  0.54  \P\     , 


(11-46) 
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(11-47) 

The  variables  cj+  and  w"  are  random  values  obtained  from  a 
normal  probability  distribution  with  zero  mean  and  standard 
deviation,  aw.  As  noted  by  Pielke  (1984,  p.  179),  this  method 
is  not  entirely  satisfactory,  since  it  depends  on  the 
particular  random  number  generator  routine. 
3.   The  NPS  Dispersion  Model 

As  mentioned  above,  L,  the  Monin-Obukhov  length  was 
proscribed  for  this  study,  while  V,  the  mean  windspeed  at  a 
height  of  two  meters  was  set  to  3m/ sec  and  the  surface 
vegetative  canopy  roughness  length  was  assumed  to  be  0.05m. 
Other  flow  parameters  required  by  the  McNider  and  NPS  particle 
models  were  computed  from  recently  developed  boundary  layer 
parameterizations,  given  here  and  by  Panofsky  and  Dutton 
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Figure  10.  Probability  density  function  of  vertical  velocity 
in  the  convective  boundary  layer.  From  Weil,  Dispersion  in 
the  Convective  Boundary  Layer,  Lectures  on  Air  Pollution 
Modelling,  Venkatram  and  Wyngaard,  1988. 
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(1984),  Sorbjan  (1990)  and  Kamada  (1992a).  More  detailed 
descriptions  of  the  underlying  theory  will  be  published  in  a 
forthcoming  article.  Many  of  the  algorithms  listed  here  are 
actually  a  part  of  the  mesoscale  windf low  simulation  which  is 
also  used  to  drive  the  McNider  model.  From  this  windf  low 
simulation  the  McNider  model  obtains  the  following:  the 
turbulent  diffusivity,  K,,,,  gradient  Richardson  number,  Ri, 
vertical  windshear,  the  windspeed  at  height  z,  vz/  and  surface 
layer  friction  velocity,  u.. 

u.  is  computed  from  L,  using  the  intermediate  Businger 
function,  \j/m,  as  determined  by  equation  11-24,  and  the  square 
root  of  the  surface  drag  coefficient, 

-i/2  _      k 


In 


1   z\  (11-48) 


so  that 

U,  =  MAX  (C^2, 0.01)   .  (11-49) 

Given  L  and  u.,   one  can  compute  the   surface  vertical 
temperature  flux, 

V¥0   =  -  ^-      .  (11-50) 

0     kgL 

This  allows  an  estimate  of  the  free  convective  boundary  layer 
turbulent  velocity  scale  (for  L  <  0) , 
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v*  =  \9Zi 


e  J 


1/3 


(11-51) 


This  can  be  modified  to  include  shear  induced  surface  layer 
turbulence  (forced  convection)  according  to 


0.4  -  — 
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z. 
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1  -  150 
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-  1 


(11-52) 


0.4 


(Kamada,  1992a) . 

For  z/L  <  -0.5,  i.e.,  above  the  surface  layer,  the 
following  forms  were  used.  The  vertical  velocity  variance  was 
given  by 


\\  =  wtt 


1.1 
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(11-53) 

where  R  =  0.2  and  D  =  0.1  are  ratios  of  the  inversion/surface 
temperature  flux  and  entrainment  zone/boundary  layer  depth 
(Sorbjan,  1990) . 

In  a  standard  atmospheric  windflow  model  with  second 
order  turbulence  closure  scheme,  the  turbulence  kinetic  energy 
(tke)  is  computed  prognostically.  If  aw2  is  supplied  as  above, 
the  horizontal  variances  follow  by  subtraction.  Without  such 
a  model,  au2  and  av2  are  obtained  diagnostically  as  follows. 

The  Andre  (1978)  third  order  turbulence  closure  used 
the  following  prognostic  form  for  vertical  velocity  variance, 
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_  2e  (11-54) 
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Assuming  steady  state  and  neglecting  diffusion,  this  can  be 
truncated  to  estimate  the  anisotropy,  crw2/e  where  e  is  the  tke. 
Mason  and  Thompson's  (1987)  large  eddy  simulation  of  neutrally 
stable  flow  showed  that  crw2/e  =  1/3  near  the  surface  and 
increased  with  height.  This  occurs  because  the  surface 
restricts  vertical  motion.  Also,  most  boundary  layer 
turbulence  results  from  "surface  no-slip"  induced  shear  which 
creates  mostly  horizontal  tke.  The  increase  with  height 
occurs  because  the  shear  drops  rapidly,  so  crw2/e  gradually 
approaches  the  isotropic  2/3  value  through  pressure  re- 
distribution, consistent  with  Mason  and  Thompson's  results. 

Like  the  tke  on  which  it  depends,  e,  the  molecular 
dissipation  rate  of  tke  also  decays  gradually  with  height.  For 
convective  boundary  layers,  Sorbjan  (1990)  parameterizes  e  as, 


3(1    -   z/zj 
<pe   =  +  Rz/Zi        ,   and  (11-55) 


e  =  0ev.//z,.   .  (11-56) 

Putting  the  above  together,  the  tke  anisotropy  ratio  can  be 
parameterized  as, 
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aj  2e0  -    e  Bw'B'z 


(11-57) 


3e0 


(Kamada,  unpublished) .   The  first  term  on  the  right  hand  side 
accounts  for  the  height  dependence  under  neutral  stability, 
while  buoyancy  flux  in  the  second  term  accounts  for  non- 
neutral  stability.   This  form  actually  corresponds  to  the 
Lenschow  et  al.  (198  0)  measurements  better  than  the  Andre 
model  itself  or  derivatives  thereof  (Therry  and  Lacarrere, 
1980) . 

So  the  horizontal  turbulent  velocity  variances  become 

°u    =  (5/9)aw2(  2e/aJ   -  1)         ,    and        (II-58a) 
a2  =   0.8au2         .  (II-58b) 

This  allows  one  to  compute  e,  the  turbulence  kinetic  energy 
(tke)  as, 

e  =  (   ou2   +  a2  +  aj  )  /2        .  (11-59) 

The  molecular  tke  dissipation  length  scale  is  parameterized  as 

1(  =   l/(l/z   +   1.4ee1/2  )  (11-60) 

from  which  the  turbulent  mixing  length  is  estimated  as, 

lk  =  MIN(31(   o2/e,    z)  .  (11-61) 
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This  finally  yields  the  required  turbulent  diffusivity, 

Km   =    0.44lkem         ,  (11-62) 

(Kamada,  1992b) .   The  buoyancy  gradient  for  the  unstable 
boundary  layer  is  given  by 

BdQ/dz   =   0  .6Q./zt(  (l-z/Zi)*3/  (z/zj1'3   - 

2R2/3(z/zi)2/3/(l    -   z/z(  +  d)2/3) 

(11-63) 

(Sorbjan,  1990) .   The  above  formulations  were  applied  to  the 

convective  (unstable)  boundary  layer  (z/L  <  -0.5)  above  the 

surface  layer. 

For  the  unstable  surface  layer  (L  <  0, 
but  z/L  >  -0.5) , 

0m  =    (1    -   28z/L)-1/4        ,  (11-64) 

Km   =  ku.z/<pm         ,  (11-65) 

e   =    (5.6KJZJ2        ,  (11-66) 

Ri    =   z/L         ,  (11-67) 

02    =    (    12    -    0.5ZJL   )m         ,  ..  (11-68) 

02   =   07  ,  (11-69) 

0i   =    (    1    -    14Z/L    )m         ,  (11-70) 

0e    =    (1    +    .75*(-z/L)2/3)3/2       ,  (H-71) 

44 


and  finally 

€0  =  uj<pjzi        ,  where  the  subscript  0  refers      (11-72) 
to  the  surface  value.  (Sorbjan,  1990) . 

For  the  neutral  to  stable  boundary  layer  and  surface 
layer,  where   L  £  0, 

a  =  2    -   10 /L        ,   and  (11-73) 

j8  =  3    -   20/L         ,  (11-74) 

(Kamada,  1992b) .   So  that, 

<Pt   =  2.2(1    -   z/zj*'2        ,  (11-75) 

02  =  0i    ,   and  (11-76) 

<p3  =   1.6(1    -   z/Zi)^2         .  (11-77) 

The  local  friction  velocity  and  Monin-Obukhov  lengths 
were  characterized  in  Sorbjan  (1990)  as 

u,  =  u.(l    -  z/z^^2        ,    and  (11-78) 

L,  =  L(l    -  z/zj3"'2-*         .  (11-79)' 

and  the  temperature  flux  at  height,  z,  was  given  by 


w'Q'z  =  w'Q'0(l    -  z/z^        .  (11-80) 

The  Brunt-Vaisala  frequency  at  height,  z,  was  given  by 


45 


BVF   =    4.3u.Q/(l    +    3.7z/L)(l    -    z/zi)p+a/z  .  (11-81) 

The   tke   dissipation  rate   at   z   was   parameterized   as 

<p(   =   3,6(1    +   3.7z/L)(l    -   z/zi)fi/z         ,    and  (11-82) 

e    =  (peuj        .  (11-83) 

The  Richardson  number  was  estimated  as 


Ri   =  ,  (11-84) 

5z   +  L 


and  the  momentum  diffusivity  was  given  roughly  by, 


ku.z(l   -  z/zt) 

Km   =  ,  (11-85) 

1+5    z/L 


(Sorbjan,  1990) .   Both  the  unstable  surface  layer  and  neutral 
to  stable  boundary  layer  cases  used: 

ou  =   u.07    ,  (II-86a) 

ov  =   au         ,  (II-86b) 

aw   =  U.0J    ,  (II-86c) 

i>h   =  21n[    (1    +    (1    -  14z/L)m)  /2  ]         ,    and        (11-87) 

0,  =  9  +  e.f  ln(z/z0)    -  \ph  )  /k  ,    where          (11-88) 

9.  =  w'Q'o/u,      .  (11-89) 
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The  mean  wind  at  height  z  and  its  shear  were  estimated  to  be 
Vz  =  u.(   ln(z/z0)    -   $m  )/k        ,  and  (11-90) 

dVz/dz   =  u.(l    +   4.7z/L)a/2(l    -   z/zj^/kz      ,  (11-91) 

(Panofsky  and  Dutton,  1984  and  Sorbjan,  1989  ) . 

The  following  details  only  the  significant  items 
distinguishing  the  NPS  and  McNider  models.   The  NPS  particle 
model  utilized  the  above  formulations  for  the  velocity 
variances  and  tke  rather  than  those  from  McNider.  For  the  NPS 
particle  model,  the  Lagrangian  vertical  timescale  was  also 
estimated  differently  according  to 

Tiw   =  U^       (  L  >   0   )  ,   and  (11-92) 

Tlw   =  0.3Zi/w.s  (   L   <    0    )  .  (11-93) 

Recognizing  that  the  horizontal/vertical  eddy  aspect 
ratio  decreases  with  increasing  stability  in  the  convective 
boundary  layer,  the  horizontal  Lagrangian  integral  timescales 
were  estimated  by 

Tlu   =    (   2    -   40 /L)    Tlw         ,       and  (11-94) 

Tlv  =   Tlu         .  .  (11-95) 

The  McNider  algorithms  for  horizontal  integral 
timescales  were  retained  for  stable  cases. 
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Rather  than  use  the  McNider  formulations,  the  vertical 
velocity  skewness  was  simulated  by  converting  every  fifth 
updraft  into  a  downdraft,  while  accelerating  the  updrafts  and 
decelerating  the  downdrafts.   So  the  updraft/downdraft 
probability  ratio  was  set  to  60%/40%/  while  the  updrafts  were 
50%  faster  than  the  downdrafts,  in  tune  with  atmospheric 
measurements . 

Non-stochastic  buoyant  forcing  was  also  added  to  the 
NPS  particle  model  for  non-zero  density  gradients.  For  pure 
buoyancy  driven  motion,  the  force  balance  is  just 

&--*"  '  (II"96) 

where  w   is  vertical  velocity,  and  69    is  the  change  in 
potential  temperature  due  to  a  displacement,  6z ,  away  from  the 
particle's  neutral  buoyancy  height  in  a  domain  where 
89/ dz   ^  0.      Here  the  neutral  buoyancy  height  is  taken  to  be 
the  initial  release  height  of  the  particle.   Then, 

dw  =  -  -f  56  dt     .  (11-97) 

u 

If  the  particle's  "memory"  were  perfect,  then  after  a 
time,  t,  its  buoyancy  forced  velocity  would  simply  be 


w 


=  -  _£  jf  50  dt     .  (11-98) 
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However,  to  simulate  dilution  by  ambient  fluid,  the  particle 
must  gradually  forget  its  neutral  buoyancy  level.   Its 
mnemonic  e-folding  time  will  be  Tlw/  the  integral  eddy 
coherence  time  scale  already  used  to  compute  the  stochastic 
component.  So  at  each  time  step  the  particle's  memory  dims  by 
the  factor, 

TT  (11-99) 

M  =  e    lw     .  K  ' 

In  this  case,  after  one  time  step, 

*i   =  -(g/e)se1st     ,  (ii-ioo) 

as  before.   But  for  subsequent  time  steps,  the  solutions  are 
w2  -  -(g/Q)St(Mse1  +  se2)       ,  (ii-ioi) 

w3  =  -(g/&)  st(M2SQ1  +  Mse2  +  se3)       ,         (11-102) 

or  in  general, 

m 
Wm   =   -(g/Q)St   E  59,.  Mm+1-1         .  (11-103) 

1 

This  series  grows  quickly  with  m.   However,  also  observe 
that 

wi+l   =  Mwt   -    (g/Q)  5ei+1St         .  (11-104) 

So  only  the  last  velocity  need  be  retained.   The  buoyant 
forcing  then  added  to  the  standard  stochastic  component 
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(adjusted  for  vertical  velocity  skewness)  to  obtain  the  total 
vertical  velocity  fluctuation. 

4 .   Random  Number  Generator  Kernel 

The  dispersion  models  utilize  normally  distributed 
random  numbers  to  generate  the  fluctuating  velocity.  These 
fluctuations  are  scaled  by  the  standard  deviation  of  the 
velocity  components  to  relate  them  to  the  turbulent  kinetic 
energy  of  the  flow  field.  The  method  used  to  generate  random 
numbers  in  this  study  was  a  variation  of  the  congruence 
method,  and  is  outlined  in  program  McNider,  subroutines  NORNG 
and  STRNUM  (Appendix  C)  .  To  check  the  distribution  of  the 
random  number  generation  routine,  the  count  distribution  was 
plotted  versus  standard  deviation  [Figure  11] .  It  appears  to 
have  a  generally  Gaussian  shape  with  perhaps  slight  skewness 
to  the  left  of  center  for  the  3600  iterations  used  in  this 
test. 
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Figure   11.    Random  number  distribution,    Program  McNider 
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III.    ANALYSIS 

A.  EQUIPMENT/ SOFTWARE 

The  logistics  difference  equation,  Henon  equations,  and 
the  McNider  Dispersion  model  were  all  written  in  standard 
FORTRAN-77,  and  run  on  Intel  386-based  personal  computers 
using  the  PROFORT  compiler.  All  graphs  were  produced  with 
Golden  Software's  Grapher  and  Surfer  programs. 

B.  LOGISTICS  DIFFERENCE  AND  HENON  EQUATIONS 

1.   Methodology  in  analyzing  the  Logistics  Difference 
Equation 

Computer  generated  data  from  the  logistics  difference 
equation  was  produced  using  Program  Feigenbm  (Appendix  1)  . 
Since  use  of  the  logistics  difference  equation  involves 
recursion,  an  iterative  series  rather  than  a  time  series  is 
created,  such  that  the  fractal  lengths  were  redefined  using 


L(i)   =  -  f   |F(i+e)  -F(i)  I  di    ,  (IH-1) 

0 

which  is  approximated  numerically  by 

N 

L(i)    =  £  \F(i   +  k)  -FU)\,  (in-2) 

where 
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F(i)    =  value  of  the  function  at  the  ith  iterate, 

i  =  iteration  step, 
and  k  =  resolution  in  terms  of  number  of  iterations. 

Analogous  to  the  Lagrangian  particle  models  studied 
later,  the  resulting  iterative  series  represents  a  single 
particle,  shifting  position,  with  each  successive  iteration 
corresponding  to  a  single  time  step.  To  insure  independence 
from  initial  values,  the  program  was  run  for  3700  iterations 
and  the  first  100  points  were  discarded.  N  was  set  to  3600 
iterations  for  all  plots,  giving  3600  point  data  sets. 

Since  recursion  does  not  permit  non-integer 
iterations,  the  values  of  k  had  to  be  integer  divisors  of 
3600.  The  45  available  values  of  k  used  are  listed  as  an  input 
data  file  in  Program  Feigenbm  (Appendix  1) . 

From  a  series  of  k  values  for  L(k)  for  a  given  ju,  a 
standard  linear  regression  routine  then  determined  DA.  The 
standard  deviation  of  DA,  was  very  small  in  chaotic 
situations,  and  large  for  during  period  doubling.  The  reasons 
for  this  are  discussed  below. 

In  computing  DA,  the  last  three  values  of  L(k)  and  k 
were  discarded.  On  most  plots,  the  curve  flattened  for 
K  >  N/3.  This  was  also  noted  by  DeCaria  (1992)  in  analyzing 
time-series  data.  So  in  the  parabolic  logistics  difference 
equation,  one  reason  for  discarding  the  highest  k  values  is 
the  rising  number  of  fold-crossings  with  larger  k.   The 
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"particle"  trajectory  is  confined  to  the  unit  interval  by  the 
left-right  symmetry  or  fold  at  mid-parabola.  Therefore,  the 
distance  traversed  during  a  fold  crossing  is  not  fully 
accounted  as  with  a  system  without  folds.  The  change  in 
"apparent"  traversed  distance  will  be  lessened  as  more  fold 
crossings  are  involved.  This  results  in  decidedly  non-linear 
slopes  for  log  L(k)  versus  log  k  for  k  which  are  a  sizeable 
fraction  of  the  window. 
2.   Analysis 

For  values  of  n  corresponding  to  a  fixed  point,  the 
plot  of  Log10L(k)  vs  Log10(k)  is  fairly  flat,  and  the 
corresponding  DA  value  is  small  [Figure  12].  For  values  of  /x 
corresponding  to  period  doubling  (period  2) ,  the  plot  shows 
large  jumps  from  a  baseline,  corresponding  to  multiple 
harmonics  of  2  [Figure  13].  The  baseline  corresponds  to  a  zero 
length.  For  these  plots,  the  zero  was  adjusted  by  adding  one 
to  avoid  a  baseline  at  -  «,  a  result  of  Log10(0)  .  This  shift  by 
one  unit  was  important  for  later  calculations;  without  it,  DA 
-*  oo  for  all  periodic  functions. 

A  period  4  plot  still  shows  regularity  [Figure  14]. 
There  appear  to  be  three  sets  of  overlapping  slope  patterns: 
one,  a  baseline  of  slope  zero;  the  other  two  with  similar 
slopes  but  different  amplitudes.  Again,  this  can  be  ascribed 
to  harmonics.  The  length  is  zero  for  every  fourth  data  point, 
and  so  will  remain  zero  for  harmonics  having  k  values  of  4 ,  8, 
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Figure  12.  Logistics  Difference  Equation,  xn+1  =  2.8xn(l-xn), 
x0  =  0.01,  Fixed  Point  Attractor. 
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Figure  13.  Logistics  Difference  Equation.  xn+1  =  3.0xn(l-xn), 
x0  =  0.01,  Period  2. 
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16,  etc.  The  length  is  nonzero  for  every  second  data  point, 
not  including  multiples  of  four.  There  is  similarity  in 
lengths  for  k  values  of  2,  6,  10,  18,  etc.  Likewise,  there  is 
similarity  in  lengths  for  k  values  which  are  not  multiples  of 
2,  e.g.,  1,  3,  5,  etc. 

At  the  onset  of  chaos  (ju  =  3.569946)  as  in  Figure  15, 
all  L(k)  lift  off  the  baseline;  without  periodicity,  all 
lengths  will  remain  non-zero.  Thus,  DA  will  increase  while  aDa 
decreases.  Also  noted  is  a  sensitivity  to  initial  conditions 
for  identical  values  of  \l ,  another  characteristic  of  chaos. 
Figure  16,  with  x0  =  0.04,  is  different  from  Figure  14  with  x0 
=  0.1.  However,  Figure  15,  with  x0  =  0.1,  is  virtually 
identical  to  Figure  17,  with  x0  =  0.01.  There  may  be  some 
intrinsic  pattern  to  this  sensitivity  to  initial  conditions. 
The  fold  crossing  patterns  might  provide  further  insight. 
However,  this  question  strays  from  the  present  focus. 

For  full  chaos,  aDa  becomes  relatively  small,  and  the 
plot  becomes  almost  linear,  as  shown  in  Figure  18.  This  plot 
highlights  why  the  last  three  points  were  discarded  in 
computing  DA;  the  plot  remains  very  linear  sans  large  values 
of  k.  The  main  reason  is  the  fold  crossing  patterns  mentioned 
earlier.  Additionally,  not  enough  data  points  are  retained  in 
the  length  computation  for  large  values  of  k  to  maintain 
similarity.  For  k=1200  and  N=3600,  the  computed  length  is  the 
sum  of  only  three  distance  measurements,  which  is  not  enough 
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Figure  14.  Logistics  Difference  Equation,  xn+1 
x0  =  0.01,  Period  4 . 
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Figure    15.    Logistics   Difference   Equation, 

xn+1    =    3 .  56994568xn(l-xn)  ,    x0   =    0.1,    Onset    of    Chaos. 
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Figure  16.  Logistics  Difference  Equation, 
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to  develop  a  good  statistical  mean  length. 

A  plot  of  DA  versus  ju  shows  many  of  the  features  noted 
in  the  previous  bifurcation  plot  [Figure  19].  At  \s.  «  3.0  DA 
jumps  from  zero  (corresponding  to  a  fixed  point)  to  ~  0.58, 
corresponding  to  period  2.  At  /i  «  3.45  is  seen  a  transition  to 
period  4.  Higher  \i  values  display  regions  of  periods  8  and  16. 
However,  at  this  level  of  resolution  period  3  2  or  higher  order 
bifurcations  cannot  be  discerned.  Within  the  chaos  regime, 
about  \x  «  3.6,  DA  seems  to  oscillate,  and  shows  several  sharp 
peaks  and  dips.  The  apparent  window  at  \l  ~  3.8  5  corresponds  to 
a  period  4  oscillation. 

A  plot  of  entropy  versus  /x  shows  similar  features 
[Figure  20].  Regions  of  period  2,  4,  and  8  are  readily 
apparent.  As  in  the  DA  versus  ju  plot  [Figure  19],  differences 
between  period  16  and  above  and  chaos  cannot  be  discerned  at 
this  level  of  resolution.  The  large  window  of  low  periodicity 
at  fi  ~  3.85  is  also  seen,  as  are  several  other  windows  of 
periodicity  at  roughly  n   «  3.63  and  ju  «  3.73. 

The  Lyapunov  exponent  also  shows  the  trends  seen  in 
the  fractal  dimension  and  entropy  [Figure  21] .  X  dips  sharply 
at  \x  «  3.25,  3.5,  and  3.55,  the  centers  for  periods  2,  4,  and 
8.  Again,  changes  beyond  period  16  are  hard  to  resolve.  In  the 
period  doubling  region,  X  =  0  indicates  that  the  function  is 
transitioning  to  the  next  bifurcation,  a  fixed  point 
fissioning  to  two  fixed  points.  For  X  >  0,  the  function  is 
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Figure  18.  Logistics  Difference  Eguation, 
xn+1  =  3.9996xn(l-xn)  ,  x0  =  0.01,  Chaos. 
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Figure  20.  Shannon  Entropy  Analysis  of  Logistics  Difference 
Equation,  xn+1  =  juxjl-xj  . 
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Figure  21.   Lyapunov  Exponent   (X)   Analysis  of  Logistics 
Difference  Equation,  xn+1  =  /xxn(l-xn)  . 


66 


chaotic.  There  are  also  numerous  "windows  of  periodicity"  in 
the  chaotic  region  which  appear  in  the  bifurcation,  fractal 
dimension,  and  entropy  plots.  These  windows  are  more  readily 
identified  by  the  Lyapunov  exponent;  if  X  drops  below  zero, 
the  function  is  periodic. 

Figure  22  shows  the  strong  correspondences  between 
entropy,  fractal  dimension,  and  Lyapunov  exponent. 
Periodicities  and  chaos  are  apparent  with  each  metric,  as  are 
some  differences.  DA  sometimes  peaks  where  S  and  X  dip,  as  at 
\i  ~  3.725.  The  low  values  of  X  and  S  indicate  periodic 
regions,  for  which  DA  sometimes  sharply  increases.  The 
bifurcation  map  [Figure  6]  provides  no  clue  as  to  the  cause. 
However,  if  the  fixed  points  are  more  widely  separated  than 
the  average  distance  between  successive  points  in  the 
surrounding  chaos  regions,  DA  will  become  relatively  large. 
This  is  less  likely  for  real  fluids  because  motion  is  then 
constrained  by  physical  forces  and  conservation  laws.  However, 
it  points  to  a  possible  schism  between  diffusion  and 
dispersion  rates  even  without  mean  flow. 
3.   Henon  Function  Analysis 

Since  the  Henon  equations  are  a  2-D  extension  of  the 
logistics  difference  equation,  similar  behavior  is  expected. 
With  two  dimensions,  the  length  measurement  is  modified  to 
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In  the  period  2  case  [Figure  7],  a  plot  of  Log10L(k)  vs 
Log10k  is  very  similar  to  the  logistics  difference  eguation  for 
period  2  [Figure  23].  Again,  for  every  k,  evenly  divisible  by 
2,  the  length  drops  to  zero.  There  are  two  sets  of  slope 
patterns:  a  baseline  with  slope  zero,  and  an  upper  non-zero 
slope  (corresponding  to  non-even  k-values) .  This  upper  slope 
appears  to  be  quite  straight,  except  for  the  largest  k  value. 

At  the  onset  of  chaos  [Figure  24],  again  the  baseline 
lifts,  and  the  overall  deviation  in  slope  is  smaller 
[Figure  25].  It  is  also  noted  that  Log10L(k)  only  approaches 
zero  when  k  is  quite  large.  As  noted  before,  for  large  k  there 
are  not  enough  distance  measurements  to  adequately  define  a 
good  statistical  mean  length;  nor  are  there  enough  points  to 
adequately  define  a  periodicity. 

At  full  chaos  [Figure  7],  the  plot  of  Log10L(k)  vs 
Log10k  is  so  linear  that,  in  regions  of  full  turbulence  or 
chaos,  DA  can  be  defined  with  only  two  data  points,  i  and  o, 
corresponding  to  one  small  inner  and  one  large  outer  scale  of 
k.  [Figure  26].  Then  DA  in  this  case  is  given  by 
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Figure  23.  Henon  Function,  a  =  0.8,  b  =  0.3,  DA  =  0.664, 
Period  Doubling. 
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This  shows  immediately  that  an  increase  in  DA  implies  a 
relative  increase  in  the  apparent  magnitude  of  the 
fluctuations  discernible  at  higher  resolutions,  i.e.,  a  shift 
in  the  amplitude  spectrum  toward  smaller  scales.  This 
interpretation  of  DA  will  be  of  use  in  analyzing  the 
atmospheric  Lagrangian  particle  models. 

There  appears  to  be  a  close  correspondence  between 
entropy  S  and  DA  for  the  Henon  function  [Figure  27].  Both  show 
jumps  corresponding  to  periods  2,  4,  8,  and  perhaps  even  for 
period  16.  Both  plots  show  corresponding  changes  when  there  is 
periodicity  within  the  chaos  region. 

A  higher  resolution  look  at  the  region, 
1.052  <  a  <  1.082,  displays  some  interesting  complexity 
[Figure  28],  Not  only  are  there  normal  bifurcations,  but  at 
a  «  1.062  a  new  branch  appears  with  no  obvious  connection  to 
previous  points.  This  does  not  correspond  to  the  fissioning 
process  observed  earlier,  but  rather  to  an  entirely  new  set  of 
solutions.  At  a  ~  1.080,  points  from  this  new  branch  break 
away  and  drift  toward  the  original  branch. 
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In  the  region,  1.052  <  a  <  1.082,  DA  jumps  in 
amplitude  whenever  the  bifurcation  map  shows  a  window  of 
periodicity  [Figure  29].  DA  based  on  x  rather  than  the  total 
distance  r  portrays  much  but  not  all  the  same  information: 
jumps  correspond  to  most  of  the  same  windows  of  periodicity 
[Figure  30].  The  strong  peaks  in  DA  correspond  to  sharp  dips 
in  S;  however,  sharp  dips  in  S  do  not  necessarily  correspond 
to  sharp  peaks  in  DA  [Figure  31]  .  So  a  jump  in  average 
distance  traversed  between  iterations  always  accompanies  a 
drop  in  position  randomness  but  not  vice  versa.  This  suggests 
that  periodicity  in  the  Henon  system  sometimes  results  from 
the  sudden  appearance  of  fixed  points  which  are  as  closely 
spaced  or  more  closely  spaced  than  the  mean  distance  between 
successive  iterations  in  the  surrounding  chaos. 

Recognizing  that  there  are  two  free  parameters  in  the 
Henon  eguation  (a  and  b)  ,  DA  and  S  were  mapped  for  values  of 
0  <  a  <  1.5  and  0  <  b  <  1.0  [Figures  32,  33].  The  two  3-d  maps 
are  remarkably  similar,  and  suggest  that  exp(DA)  would  be  of 
the  same  order  as  S.  The  regimes  of  low  periodicity  are  well 
defined.  Where  DA  and  S  are  zero,  the  Henon  function  has 
period  1,  corresponding  to  one  fixed  attractor.  The  first 
jumps  to  periods  2  and  4  are  well  defined  at  this  level  of 
resolution.  There  are  other  steps  of  discernible  width  for 
high  a  and  b  values,  on  the  left  and  back  sides  of  the 
graphical  "mountain".  Another  interesting  feature  is  the 
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Figure  29.   Fractal   Dimension   (DA) 
Function,  testing  Ar  with  b  =  0.3. 
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Figure  30.   Fractal   Dimension   (DA)   Analysis   of   Henon 
Function,  testing  Ax,  with  b  =  0.3. 
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behavior  for  high  values  of  b,  where  both  DA  and  S  show  sudden 
increases  in  value. 

4.   Conclusions  drawn  from  Logistics  Difference  and  Henon 
equation  analysis 

The  self-affine  fractal  dimension,  DA,  is  a  good 
discriminator  of  low  frequencies  in  data,  e.g.,  periods  2,  4, 
and  8  of  a  given  length  data  set.  However,  at  high 
frequencies,  period  16  and  above  in  the  initial  bifurcation, 
it  is  difficult  to  use  DA  to  differentiate  between  periodicity 
and  turbulence  or  chaos.  The  entropy,  S,  seems  related  to  DA, 
but  S  measures  the  evenness  of  particle  position  distribution, 
while  DA  measures  the  changing  apparent  jaggedness  of  the 
function  as  the  resolution  varies.  A  sharp  change  in  S  always 
accompanies  a  sharp  change  in  DA,  but  not  always  vice  versa. 
Therefore,  DA  and  S  generally  show  the  same  trends  within  the 
chaos  region,  but  not  always.  In  such  cases  diffusion  and 
dispersion  rates  are  not  equivalent. 

The  Lyapunov  exponent  provides  perhaps  the  most  definitive 
information:  for  a  given  value  of  X,  we  know  for  certain 
whether  the  function  is  stable,  periodic,  or  chaotic. 
Additionally,  the  Lyapunov  exponent  should  be  able  to  describe 
the  degree  of  chaos:  the  greater  the  value  of  X,  the  faster 
the  orbital  trajectory  diverges,  and  therefore  the  more 
chaotic  the  function.  Unfortunately,  the  Lyapunov  exponent  is 
not  readily  determined  for  multi-dimensional  systems. 
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Figure   32.       Fractal       Dimension       (DA) 
Function,    xn+1   =    l-axn2+yn,    yn+1   =   bxn. 


Analysis       of      Henon 


82 


<0 


Figure   33.    Entropy    (S)    Analysis   of   Henon   Function, 
xn+i  =   l~axn2-yn,    yn+1  =  bxn. 
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Analysis  of  these  two  simple  systems,  the  logistics 
difference  and  Henon  sets,  shows  that  S,  DA,  and  X  display 
correlated  but  not  identical  behavior,  since  they  measure 
rather  different  things.  These  insights  can  now  be  applied  in 
studying  atmospheric  Lagrangian  particle  models. 

C.   MCNIDER  PARTICLE  DISPERSION  MODEL  ANALYSIS 
1.   Methodology 

Four  aspects  of  the  McNider  Particle  Dispersion  Model 
were  studied:  the  divergences  of  vertical  position,  total 
position,  velocity,  and  phase  velocity.  The  position  was 
measured  in  radial  distance  from  the  particle  origin,  which 
was  arbitrarily  set  at  x,y  =  0,  z  =  100  meters.  To  simulate 
the  real  atmosphere,  the  boundary  layer  inversion  height  was 
varied  linearly  from  2,000  to  200  meters  as  1/L,  the  inverse 
Monin-Obukhov  length,  was  varied  from  -0.2  to  0.1.  Mean 
velocities  were  set  to  zero  and  total  particle  velocity 
reflection  was  assumed  at  the  ground  surface  and  inversion. 

Hints  as  to  model  behavior  are  again  provided  by  the 
bifurcation  maps,  which  now  depict  the  expansion  and 
distribution  of  particle  range  with  decreasing  stability,  as 
measured  by  1/L.  Model  performance  was  checked  against 
standard  geophysical  measures  such  as  the  Brunt  Vaisala 
Frequency  (BVF) ,  turbulent  kinetic  energy  (tke) ,  and  vertical 
velocity  variance  (aw2)  ,  as  well  as  against  the  two  readily 
calculated  chaos  measures.,  entropy,  S,  and  the  self-affine 
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fractal  dimension,  DA,  for  the  range  of  1/L  values  seen  in  the 
atmosphere.  Using  100  partitions,  the  computed  maximum 
possible  entropy  for  all  plots  in  the  analysis  was 

S^-4.6   .  (III-5) 

For  this  study  1/L  was  set  at  a  particular  value,  a 
single  particle  was  released  and  followed  for  3,600  time 
increments  of  1/6  seconds  each  (600  seconds  total) ,  then  1/L 
was  incremented  and  the  analysis  repeated,  until  the  entire 
range  was  spanned  using  1/L  increments  of  0.01. 
2.   McNider  No-Skew  Routine 

For  real  unstable  atmospheres,  the  vertical  velocity 
distribution  is  negatively  skewed  [Figure  10] .  Updrafts  have 
higher  velocities  and  thus  occupy  less  volume  than  downdrafts. 
However,  for  the  initial  analysis  of  the  McNider  model,  the 
vertical  velocity  turbulence  distribution  was  not  skewed.  This 
was  done  to  check  the  model  without  the  ad  hoc  method  McNider 
developed  to  introduce  skewness,  and  which  as  outlined  in 
Pielke  (1984,  pp.  178-179)  appears  to  be  incorrect. 

An  r* versus  1/L  plot  of  the  particle  position  shows 
that  on  the  negative  (unstable)  side,  DA  and  S  follow  roughly 
the  same  trend  up  to  1/L  =  0  [Figure  34].  From  the  analysis  of 
DA  in  the  Henon  system,  this  suggests  that  the  ratio  of  large 
to  small-scale  vertical  distances  between  points  varies  little 
in  the  unstable  regime  without  overlaying  skewness.  However, 
DA  again  jumps  suddenly  when  1/L  exceeds  zero.  Since  vertical 
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fluctuations  tend  to  decay  in  general  with  increasing 
stability,  one  expects  that  S  would  decay  and  DA  would  tend 
upward  steadily  from  left  to  right  for  real  turbulence  in  the 
atmosphere.  Again,  the  discontinuities  near  1/L  =  0  are 
probably  due  to  discontinuities  across  the  neutral  transition 
in  the  McNider  algorithms.  The  small  scale  fluctuations  in  S 
with  1/L  are  probably  due  to  the  random  nature  of  the  particle 
paths. 

On  the  positive  (stable)  side  of  1/L,  both  DA  and  S 
show  sudden  dramatic  step  increases,  followed  by  DA  tending 
upward  while  S  tends  downward.  Again,  this  highlights  the  fact 
that  S  and  DA  are  not  measuring  the  same  thing.  The  sudden 
jumps  in  S  in  the  McNider  model  seem  antithetical  to  the  fact 
that  negative  buoyancy  tends  to  suppress  vertical  motion  and 
render  it  oscillatory  in  real  fluids.  For  growing  positive 
values  of  1/L,  the  decrease  in  S  shows  that  the  distribution 
of  points  within  the  domain  becomes  less  uniform,  while  the 
rise  in  DA  suggests  that  the  distance  between  successive 
points  decreases  more  slowly  on  the  small  scale  than  for  the 
larger,  longer  time  scale  fluctuations.  This  is  qualitatively 
consistent  with  measurements  under  increasingly  stable 
conditions  in  real  fluids  (Stull,  1988) . 

The  bifurcation  diagram  [Figure  35]  also  shows  that 
for  1/L  <  0  the  total  range  of  positions  steadily  decreases 
with  increasing  stability,  while  for  1/L  >  0  it  hardly  varies. 
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This  is  because  the  total  position  change  includes  both 
vertical  and  horizontal  fluctuations,  and  horizontal 
fluctuations  are  much  less  dependent  on  stability. 

These  plots  reveal  other  weaknesses  in  the  McNider 
model.  Unrealistically  sharp  fluctuations  seem  to  occur  near 
the  neutral  stability  transition  (1/L  =  0)  ,  probably  due  to 
the  dichotomous  nature  of  the  McNider  algorithms,  one  set  for 
stable  conditions,  another  for  unstable.  From  equation  11-18, 
it  is  also  clear  for  isotropic  turbulence  that  aw2  =  2/3  tke. 
Stability  suppresses  vertical  turbulence,  while  instability 
enhances  it.  So  for  1/L  <  0,  aw2  should  exceed  2/3  tke.  Figure 
34  shows  increasingly  qualitatively  incorrect  ratios  with 
increasing  instability.  aw2  also  drops  to  near  zero  at  neutral 
stability  (1/L  =  0)  due  to  very  low  mean  windshear  values  in 
the  mesoscale  flow  simulation.  This  may  not  be  as  significant 
a  problem  in  the  prognostic  windflow  models  for  which  the 
McNider  model  was  originally  designed  as  it  is  in  the 
similarity-based  flow  model  simulator  used  in  this  study. 
However,  an  artificial  "dip"  would  probably  still  appear. 

The  bifurcation  map  of  the  total  3-D  position  seems  to 
show  a  fairly  constant  range  for  1/L  >  0,  a  jump  at  zero,  and 
a  varying  range  for  1/L  <  0  [Figure  35].  The  plot  density  is 
too  heavy  to  discern  actual  distribution  patterns  for  given 
values  of  1/L,  although  the  computer  screen  displays  a  Moire- 
like pattern.  What  can  be  extracted  from  this  plot  is  the 
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Figure  34.  McNider  Particle  Dispersion  Model.  Test  of  total 
3-D  distance,  Ar,  with  no  skew. 
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general  density  of  points  and  the  upper  and  lower  boundaries 
(in  magnitude)  of  position. 

In  contrast  to  Figure  34,  a  plot  of  only  the  vertical 
position  change  (Az)  reveals  that  the  pattern  for  DA  is  quite 
similar  to  that  for  S  [Figure  36].  The  bifurcation  map  shows 
an  opposite  jump  in  the  vertical  position  as  1/L  crosses  the 
zero-point  [Figure  37]. 

The  velocity  distribution  for  various  values  of  1/L 
seems  fairly  constant,  so  that  the  S  curve  has  only  a  slight 
upward  trend  [Figure  38],  except  at  values  near  1/L  =  0. 
However,  DA  tends  downward  strongly  from  left  to  right,  jumps 
at  zero,  then  again  tends  downward  at  around  1/L  =  0.7,  where 
it  begins  increasing  in  value.  This  upward  trend  in  DA  is 
probably  caused  by  increasing  negative  buoyancy.  The 
bifurcation  map  [Figure  39]  shows  that  the  velocity  range 
tends  generally  downward,  with  a  jump  near  1/L  =  0,  and  a 
change  in  trend  at  1/L  =  0.7. 

The  phase  velocity  plots  are  almost  identical  to  the 
total  position  plots.  This  is  because  in  calculation  phase 
velocity  the  velocities  are  small  in  magnitude  compared  to  the 
position  values. 

3.   McNider  Model,  Skew  On 

A  glance  at  Figure  42  shows  a  big  problem  with  the 
vertical  velocity  skewness  algorithm.  The  tke  and  aj  values 
are  much  too  large.  An  expected  maximum  aw2  will  be  ~  10m2/s2, 
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Figure  36.   McNider  Particle  Dispersion  Model 
vertical  position,  Az ,  with  no  skew. 
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Figure  38.  McNider  Particle  Dispersion  Model,  Test  of  total 
velocity,  no  skew. 
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Figure  40.  McNider  Particle  Dispersion  Model.  Test  of  Phase 
Velocity,  with  no  skew. 
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Figure  41.  McNider  Particle  Dispersion  Model  Bifurcation 
Map.  Test  of  phase  velocity,  with  no  skew. 
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while  the  plot  shows  ~  100  m2/s2.  It  seems  that  the  algorithm 
is  incorrect  for  very  unstable  values  of  1/L.  However,  the 
stable  side  (positive  values  of  1/L)  shows  well  behaved, 
reasonable  values  for  tke  and  aw2. 

Examining  only  the  stable  side  (positive  1/L)  of 
Figures  42-49,  DA  and  S  are  surprisingly  flat  for  velocity, 
but  not  for  total  3-D  position,  vertical  position  change,  or 
phase  velocity.  Also,  the  bifurcation  plot  of  velocity  range 
becomes  guite  small  for  positive  1/L  (Figure  47)  .  Oddly 
enough,  the  small  values  begin  not  at  1/L  >  0,  but  1/L  >  0.02. 
This  corresponds  to  a  sharp  drop  in  entropy  and  fractal 
dimension  near  this  value  of  1/L,  suggesting  that  the  skewness 
algorithm  is  responsible.  However,  re-examining  Figures  3  6  and 
38  with  "no-skew",  entropy  and  DA  drop  sharply  at  1/L  =  0.02 
as  well  as  at  1/L  =  0,  suggesting  that  something  other  than 
the  skewness  algorithm  is  responsible,  perhaps  the  inherent 
dichotomy  in  eguations  11-44  or  II-37a. 

Again,  this  sharp  transition  near  neutral  values  of 
1/L  does  not  seem  to  accurately  reflect  nature,  but  rather 
appears  to  be  a  discontinuity  in  the  algorithm  solutions. 

Figure  47  also  shows  that  the  velocity  range  is  much 
too  high  for  the  unstable  case;  an  expected  velocity  value  is 
on  the  order  of  10  m/s  or  less.  Oddly,  there  are  several 
windows  where  the  velocity  range  drops  to  low  values,  and  are 
roughly  what  is  expected,  i.e.,  for  1/L  ~   -0.475.  However,  at 
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Figure  42.  McNider  Particle  Dispersion  Model.  Test  of  total 
position,  Ar,  with  skewness. 
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Map.  Test  of  total  position,  Ar,  with  skewness. 
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Figure  44.   McNider  Particle  Dispersion  Model.   Test  of 
vertical  position,  Az ,  with  skewness. 
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Figure    45.     McNider    Particle    Dispersion    Model     Bifurcation 
Map.    Test   of   vertical   position,    Az ,    with   skewness. 
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Figure  46.  McNider  Particle  Dispersion  Model.  Test  of  total 
velocity,  with  skewness. 
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Figure    47.     McNider    Particle    Dispersion    Model     Bifurcation 
Map.    Test   of   total   velocity,    with   skewness. 
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Figure  48.  McNider  Particle  Dispersion  Model.  Test  of  phase 
velocity,  with  skewness. 
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Figure  49.  McNider  Particle  Dispersion  Model  Bifurcation 
Map.  Test  of  phase  velocity,  with  skewness. 
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more  negative  values  of  1/L  the  aw2/tke  ratio  indicates  that 
nearly  all  the  turbulence  is  in  the  vertical  component  which 
is  not  seen  in  the  atmosphere  nor  is  physically  likely  (Stull, 
1988) . 

4.   NPS  Model,  Skew  On 

The  NPS  model  employed  the  double  Gaussian  skewness 
scheme  as  described  in  the  theory  section.  It  avoids  the 
problems  seen  with  the  McNider  model  skewness  scheme  regarding 
the  inappropriate  high  values  for  velocity,  tke,  and  aw2  for 
unstable  values  of  1/L.  Also,  the  tke  and  aw2  better  reflect 
the  behavior  of  real  fluids:  aw2/tke  exceeds  2/3  for  unstable 
values  of  1/L,  but  not  for  stable  values.  The  sharp 
discontinuity  in  transition  from  unstable  to  stable  1/L  is 
reduced  somewhat. 

Curiously,  there  is  a  distinct  pattern  to  the  strong 
fluctuations  in  DA,  S,  tke,  and  ctw2  for  unstable  1/L  but  not 
for  stable  1/L  values.  This  was  also  observed  in  the  McNider 
skewness  routine.  The  cause  has  not  been  determined;  it  could 
be  due  to  limitations  in  the  random  number  generator  or 
perhaps  an  undiscovered  program  error.  This  is  most  obvious  in 
the  bifurcation  map  of  vertical  position  range,  Az  [Figure 
53],  where  there  are  sharp  jumps. 

The  phase  velocity  plots  are  similar  to  the  total 
position  plots  in  all  cases  examined.  The  reason  is  readily 
determined:  the  phase  velocity  is  dominated  by  the  total 
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position,  which  varies  from  near  zero  to  1000  meters,  while 
the  velocity  varies  from  near  zero  to  roughly  3  meters/ second. 
As  in  the  McNider  plots  there  is  little  indication 
that  the  entropy  is  affected  by  a  rising  BVF,  while  DA  seems 
to  rise  as  BVF  rises,  most  notably  in  the  vertical  position 
plots.  The  BVFs  indicate  a  long  period  compared  with  the 
1/6  second  timesteps.  The  entropy  analysis  will  not  detect  a 
particle  riding  on  a  low-frequency  wave  with  long  period.  In 
such  a  case  the  particle  distribution  over  time  may  appear  to 
be  roughly  uniform  between  a  minimum  and  maximum  value.  This 
constraint  is  similar  to  the  resolution  requirements  for 
entropy  computation  noted  earlier.  If  the  At  in  the  model  were 
increased  from  1/6  second  to  say  10  seconds,  while  still 
retaining  3  600  iterations,  the  entropy  might  also  show  a 
variance  with  the  BVF. 
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Figure  50.  NPS  Particle  Dispersion  Model.  Test  of  total 
position,  Ar,  with  skewness. 
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Figure  51.  NPS  Particle  Dispersion  Model  Bifurcation  Map 
Test  of  total  position,  Ar,  with  skewness. 
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Figure  52.  NPS  Particle  Dispersion  Model.  Test  of  vertical 
position,  Az ,  with  skewness. 
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Figure  53.  NPS  Particle  Dispersion  Model  Dispersion  Map 
Test  of  vertical  position,  Az ,  with  skewness. 
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Figure  54.  NPS  Particle  Dispersion  Model.  Test  of  velocity, 
with  skewness. 
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Figure  55.  NPS  Particle  Dispersion  Model  Bifurcation  Map 
Test  of  velocity,  with  skewness. 
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Figure  56.   NPS  Particle  Dispersion  Model 
velocity,  with  skewness. 
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Figure   57.     NPS    Particle    Dispersion    Model    Bifurcation    Map 
Test   of   phase  velocity,    with  skewness. 
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IV.  CONCLUSIONS 

Chaos  metrics  such  as  the  self -af fine  fractal  dimension, 
DA,  entropy,  S,  and  Lyapunov  exponent,  X,  are  useful  tools  in 
the  analysis  of  system  characteristics  and  behavior.  DA 
measures  the  jaggedness  of  a  trajectory,  or  relative  magnitude 
of  small  versus  large  scale  fluctuations.  The  entropy  measures 
the  randomness  of  state  distribution.  The  Lyapunov  exponent, 
difficult  to  implement  in  higher  dimensions,  measures  the 
divergence  of  trajectories,  and  leads  to  a  predictability  time 
scale.  As  in  both  the  chaos  and  atmospheric  systems,  DA  and  S 
often  correspond  closely,  but  not  always.  The  Henon  equations 
also  demonstrated  that  exp(DA)  is  often  of  the  same  order  as 
S. 

These  chaos  metrics  when  used  in  conjunction  with  standard 
geophysical  measures  such  as  turbulent  kinetic  energy,  tke, 
vertical  velocity  variance,  aw2,  and  the  Brunt-Vaisala 
Frequency,  BVF,  can  be  applied  usefully  to  atmospheric  3-D 
particle  dispersion  models.  The  McNider  model  as  listed  in 
Pielke  (1984)  has  an  inherent  weakness  in  the  skewness  of 
vertical  velocity  variance.  This  weakness  leads  to  obviously 
incorrect  values  of  tke,  aw2,  positions,  and  velocities  for 
unstable  values  of  the  inverse  Monin-Obukhov  length,  1/L.  The 
McNider  algorithms  also  display  an  inherent  discontinuity  as 
the  inverse  Monin-Obukhov  length  crosses  the  zero  point, 
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reflecting  a  trans it ion 'from  unstable  to  stable  buoyancy.  This 
sharp  discontinuity  does  not  correspond  to  the  real 
atmosphere,  which  suggests  that  model  performance  suffers  at 
small  values  of  1/L. 

The  NPS  particle  dispersion  model's  measured  performance 
seems  to  better  reflect  reality,  it  does  not  share  the  McNider 
model's  weakness  in  the  skewness  of  vertical  velocity 
variance,  and  displays  more  realistic  ratios  of  aw2  to  tke.  The 
NPS  Particle  Dispersion  Model  also  reduces  the  discontinuity 
in  transition  from  negative  to  positive  values  of  1/L. 
However,  the  NPS  model  still  has  a  problem  with  unstable 
inverse  Monin-Obukhov  lengths  (1/L  <  0)  ;  it  displays  a  pattern 
to  the  vertical  position  distribution  not  reflective  of  real 
fluid  behavior.  Since  the  NPS  model  was  developed  in  response 
to  the  current  study,  there  may  still  be  problems  in  coding  or 
in  the  random  number  generator.  The  NPS  model  requires  further 
refinement  to  model  particle  behavior  in  a  fluid 
realistically. 

Previous  performance  metrics  for  atmospheric  particle 
models  were  only  designed  to  measure  aggregate  particle 
behavior.  The  above  chaos  metrics  parameters  also  offer  some 
insight  into  the  model  behavior  of  individual  particles. 
Results  indicate  that  atmospheric  dispersion  models  require 
further  development  at  a  fundamental  level  in  replicating 
turbulent  diffusion.  For  example,  large  fractal  dimensions  in 
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atmospheric  Lagrangian  particle  mode»ls  seem  indicative  of 
strongly  stable,  rather  than  unstable,  conditions,  contrary  to 
what  seems  likely  for  real  turbulence.  This  points  to  the 
neglect  in  such  models  of  motions  at  scales  other  than  the 
dominant  scale  determined  by  <VV,W«  Thus,  atmospheric 
Lagrangian  particle  models  may  simulate  only  single  scale  or 
at  most  a  limited  range  of  large  scale  diffusion  processes. 
Fluctuations  more  in  accord  with  real  velocity  spectra  might 
display  more  realistic  entropy  and  fractal  behavior.  As  is,  S 
and  DA  behave  oppositely  at  times,  which  indicates  that 
diffusion  and  dispersion  rates  do  not  have  equivalent  meaning, 
even  in  the  absence  of  mean  windflow. 

Simple  periodic  behavior  preliminary  to  chaos  in  classical 
bifurcatory  systems  is  also  apparently  not  entirely  equivalent 
to  laminar  wave  behavior  prior  to  the  onset  of  turbulence.  For 
example,  if  the  time  resolution  of  the  analysis  is  much 
shorter  than  the  period,  the  Shannon  entropy  for  an 
atmospheric  Lagrangian  particle  model  may  be  quite  large  for 
wave  motion,  but  small  for  periodic  behavior.  This  is  because 
the  number  of  possible  states  in  wave  motion  is  not  limited  to 
the  amplitude  extrema  as  in  the  periodic  motion. 

One  weakness  of  this  study  of  Lagrangian  particle  model 
performance  is  the  unavailability  of  real  data  for  comparison 
purposes.  As  stated  earlier,  atmospheric  data  is  normally 
obtained  in  the  Eulerian  rather  than  Lagrangian  frame.  This 
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suggests  a  need  for  development  of  Lagrangian-based  sensing  in 
fluid  turbulence  experiments.  Position  and  velocity 
measurements  of  particle  markers  in  a  wide  range  of 
atmospheric  conditions,  both  stable  and  unstable,  as  a 
function  of  time,  with  rapid  response  remote  sensors,  would  be 
useful  in  developing  better  models  of  particle  dispersion  and 
diffusion.  As  yet,  since  the  gathering  of  data  from  real 
fluids  in  the  Lagrangian  frame  is  not  feasible,  these  study 
results  suggest  that  more  extended  application  of  chaos 
metrics  to  Eulerian  measurements  of  turbulence  in  real  fluids 
is  warranted  in  order  to  gauge  the  performance  of  Eulerian 
turbulence  models.  Such  applications  need  not  be  restricted  to 
small-scale  phenomena,  since  mesoscale  Rayleigh-Benard 
convection  also  displays  transitions  from  laminar  cellular  to 
chaotic  behavior  (Agee  et  al . ,  1973) .  Model  improvements  based 
on  such  tests  may  then  be  extended  to  the  Lagrangian  frame. 

With  regard  to  the  Eulerian  frame,  both  fractal  dimension 
and  Shannon  entropy  are  simple  calculations  which  can  be 
performed  real-time  in  situ  with  current  sampling  devices,  and 
should  be  simpler  to  implement  than  FFTs,  since  there  are  no 
transform  matrices.  One  area  to  consider  when  conducting 
measurements  of  the  self-affine  fractal  dimension  is  the  size 
of  both  e,  the  time  increment,  and  T,  the  time  window  width 
over  which  the  length  is  being  defined,  e  has  a  minimum  size 
dictated  by  the  response  time  of  the  sensors,  while  the  size 
of  T  is  critical  when  looking  for  intermittent  or  low 
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frequency  events,  e.g.,  gravity  waves  governed  by  the  Brunt- 
Vaisala  Frequency.  If  emin  is  too  small,  then  the  time 
increment  will  fall  in  the  range  where  the  velocity 
autocorrelation  decay  is  not  exponential  but  constrained  to  be 
parabolic  due  to  continuity  and  the  viscosity  of  real  fluids. 
At  the  same  time  e,^  must  be  at  least  -  3  orders  of  magnitude 
smaller  than  T  to  establish  statistical  validity.  Establishing 
the  appropriate  e^  is  crucial  to  distinguish  waves  from 
turbulence,  and  probably  applies  equally  for  distinguishing 
intermittent  and  coherent  structures  in  general.  (Kamada  and 
DeCaria,  1992) 

Another  area  for  further  study  is  using  the  Lyapunov 
exponent  in  analyzing  3-D  turbulence.  Although  difficult  to 
compute  in  3-D,  X  provides  a  definite  test  of  chaos  in  the 
particle  diffusion  rate.  A  study  of  all  three  chaos  metrics 
seems  necessary  to  clarify  the  relationship  of  diffusion  rate 
to  dispersion  rate  in  these  models. 

Standard  geophysical  turbulence  measures  such  as  tke  and 
aw2  might  also  be  applied  to  classical  chaotic  systems  such  as 
the  logistics  difference  and  Henon  equations.  The  definitions 
of  timescale  and  averaging  time  must  first  be  established ■, 
e.g.,  one  iteration  equals  one  time  unit.  The  corresponding 
tke  would  be  zero  for  a  fixed  point,  and  definite  magnitude 
changes  would  occur  as  the  function  bifurcates  to  periods  2, 
4,  8,  and  to  chaos. 
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The  utility  of  chaos  metrics  in  analyzing  the  3-D 
Monte  Carlo  based  particle  dispersion  models  also  suggests 
possible  utility  for  other  laminar/turbulent  phenomena.  These 
might  include  plasmas,  acoustics,  and  laser  cavities.  Solid 
state  free  electron  gas  systems  may  also  find  these  chaos 
metrics  useful  in  describing  phonon  and  electron  transport. 
Phonon  resonance  in  crystal  lattices  also  reminds  one  that 
particle  behavior  in  lattice  gases  and  cellular  automata  may 
be  studied  with  such  chaos  metrics.  These  measures  might  also 
be  useful  in  modeling  gas  or  liguid  phase  complex  chemical 
kinetics  or  radiation,  which  have  long  been  simulated  by  Monte 
Carlo  schemes. 
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APPENDIX  A.  LOGISTICS  DIFFERENCE  EQUATION  ANALYSIS  PROGRAM 

c************************************************************************ 

PROGRAM  Feigenbm 
c 

c   This  program  computes  self-affine  fractal  dimension,  Da,  for  the 
c  Feigenbaum  recursion  relation, 
c 

c       x(i+l)  =  4  lambda  x(i)(  1  -  x(i), 
c 
c 

c   as  a  function  of  lambda.   Da  is  defined  for  the  iteration  series  as: 
c 
c  d  Ln[  L(i)  ] 

c       Da  =  -  ,  where  L(i)  =  (1/i) 

c  d  Ln[  i  ] 

c 

c  N 

c  =   SUM  |F(i)  -  F(i-l) | 

c  i=k,k 

c 

c   and  k  is  an  exact  integer  divisor  of  N.   Here  we  choose 

c   N  =  3,600  to  give  us  at  least  three  decades  to  use  in  computing  Da. 

c   This  also  allows  48  exact  integer  divisors  of  N  in  the  k  array.   So 

c  this  gives  us  48  points  of  i  vs.  L(i)  in  the  linear  regression  for  Da. 

Q**************************  ********************************************** 

c 

DOUBLE  PRECISION  lambda,  lglngth(48,  100),  lgk(48),  y,  sum,  x 
DIMENSION  y (0:3600) ,  sum(48,  100),  k(48),  nyval(48),  xy(2,48) 
&  ans(16),  Dalmbda(3,100) 

DATA  k/1,  2,  3,  4,  5,  6,  8,  9,  10,  12,  15,  16,  18,  20,  24,  25, 
&      30,  36,  40,  45,  50,  60,  72,  75,  80,  90,  100,  103,  106,  109, 
&     116,  120,  144,  150,  180,  200,  225,  240,  300,  360,  400,  450, 
&  600,  720,  900,  1200,  1800,  3600/ 

OPEN  (1,  FILE=' feigenis.dat' ,  RECL=255,  STATUS= ' new' ,  ERR=700) 
OPEN  (2,  FILE=" feigenll.dat' ,  RECL=255,  STATUS= ' new ' ,  ERR=800) 
c 

y(0)  =  0.0 
c  Get  Da  for  different  values  of  lambda  from  0  to  1. 
m  -  100 

DO  500  1  =  1,  m 
el  =  1 

lambda  =  0.01*el 

IF  (lambda  .gt.  0.9999)  lambda  =  0.9999 
x       =  0.01 
c  Create  the  iteration  series  for  a  given  lambda 
N      =  7200 

WRITE  (1,*)'    i       y(i) 
DO  200  i  =  1,  N 

y(i)  =  4.*lambda*x*(l-x) 
x    =  y(i) 
WRITE  (1,150)  i,  y(i) 
150  FORMAT  (lx,  14,  3x,  f8.5) 

200  CONTINUE 

c  With  iteration  series,  get  lengths  and  self-affine  fractal  dimension 
WRITE  (2,  *)'    Ln(e)    Ln  [L(e)]' 
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c  Get  total  length  as  sum  of  differences  for  different  spacings  of  i 
DO  400  j  =  1,  48 

sum(j,l)  =  0. 

DO  300  i  =  3600  +  k(j),  7200,  k(j) 

sum(j,l)  =  sum(j,l)  +  abs(  y(i)  -  y(  i-k(j)  )  ) 
d  IF  (j  .gt.  44) 

d     &  WRITE(6,250)  y ( i) ,y ( i-k( j ) ) , sum( j , 1 ) , i, j , k( j ) 

d250  FORMAT  (  lx,  3f8.4,3l5) 

300  CONTINUE 

d  IF  (j  .gt.  44) 

d     &      WRITE(6,*) 'sum( ' , j, ' ,1)  =  ',sum(j,l) 

lglngth(j,l)  =  alogl0(  sum(j,l)  ) 

lgk(j)  =  alogl0(  float(k(j))  ) 

xy(l,j)  =  lgk(j) 

xy(2,j)  =  lglngth(j,l) 

WRITE  (2,350)  lgk(j),  lglngth( j , 1) 

FORMAT  (lx,  f8.5,  3x,  f8.5) 

CONTINUE 


350 
400 
c 

c 
c 

c 
c 
c 
c 
c 
c 


Do  linear  regression  to  get  Da  from  loglO(k)  vs.  loglO  [L(e)]  values 

CALL  STLNRG(n,  nyval,  xy,  ans) 

-ans(2)  contains  Da 

ans (14)  contains  the  standard  deviation  of  Da 

Put  lambda,  Da,  and  aDa  in  a  3  by  1  array 


Dalmbda(l,l)  =  -ans (2) 
Dalmbda(2,l)  =  -ans (14) 
Dalmbda(3,l)  =  lambda 
500       CONTINUE 

WRITE  (3,550)  ( (Dalmbda( i, 1 ) ,  i  =  1,3),  1  =  l,m) 
550    FORMAT  (lx,  3(f9.5,2x)) 

STOP 
700   WRITE  (6,*) 'error  opening  feigenis.dat' 

STOP 
800   WRITE  (6,*) 'error  opening  feigenll.dat' 
STOP 
END 
c*************************************************^ 

subroutine  stlnrg  (k,  nyval,  xy,  ans) 


This  subroutine  computes  the  slope  and  other  statistics  for  a 
linear  regression  with  several  y  values  for  each  x  value  or  with 
one  independent  variable. 

argument   use     description 

k       input   Number  of  different  x  values 

nyval    input   Array  of  number  of  y  values  for  each  x 

value.  Dimensioned  k 
xy      input    array  of  x  and  y  pairs  (xl,  yl,  x2 ,  y2,  etc.) 

of  dimension  (2,k) 

ans     output  Array  of  results  computed,  dimensioned  16. 

ans 

1  Number 
*2  Slope 
3  Mean  of  y 
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c 

4 

Mean  of  x 

c 

5 

y-intercept 

c 

6 

Sum  of  y**2 

c 

7 

Sum  of  squares  mean 

c 

8 

Sum  of  squares  slope 

c 

9 

Residual 

c 

10 

standard  deviation  of  x 

c 

11 

standard  deviation  of  y 

c 

12 

standard  deviation  error 

c 

13 

standard  deviation  y-bar 

c 

*14 

standard  deviation  slope 

c 

15 

standard  deviation  y-intercept 

c 
c 
c* 

16 

F-Ratio  for  slope 

****************************************************************** 

DIMENSION  nyval(l), 

xy(l),  ans(16) 

sumx  =  0.0 

sumy  =  0.0 

. 

sumx2  =  0.0 

ans(6)  =  0.0 

sumxy  =  0.0 

n  =  0 

nx  =  1 

ny  =  2 

DO  20  j  =  l,k 

nyval(j)  =  1 

m  =  nyval ( j ) 

n  =  n  +  m 

DO  10  i  =  1,  m 

sumx  =  sumx  +  xy(nx) 

sumx2  ■  sumx2  +  xy (nx) *xy (nx) 

sumy  =  sumy  +  xy(ny) 

ans(6)  =  ans(6)  +  xy (ny) *xy (ny) 

sumxy  =  sumxy  +  xy (nx) *xy (ny) 
10    ny  =  ny  +  1 

nx  =  nx  +  nyval ( j )  +1 
20    ny  =  ny  +  1 

en  =  n 

ans(l)  =  en 

si  =  ans(l)*sumx2  -  sumx*sumx 

s2  =  ans(l)*sumxy  -  sumx*sumy 

enl  =  en  -  1.0 

en2  =  enl  -  1.0 

ans(2)  =  s2/sl 

ans(3)  =  sumy/en 

ans(4)  =  sumx/en 

ans(5)  =  ans(3)  -  ans (2 ) *ans (4) 

ans ( 7 )  =  ans ( 3 ) *  sumy 

ans(8)  =  ans(2)*s2/en 

ans (9)  =  ans(6)  -  ans(7)  -  ans (8) 

s4  =  ans(9)/(en  -  2.0) 

endsl  =  en/sl 

ans(10)  =  sqrt(1.0/(enl*endsl) ) 

ans (11)  =  sqrt((ans(6)  -  ans(7))/enl) 

ans (12)  =  sqrt(s4) 

al3  =  s4/en 

ans (13)  =  sqrt(al3) 

al4  =  s4*endsl 

ans (14)  =  sqrt(al4) 

ans(15)  =  sqrt(al3  +  al4*ans (4) *ans (4) 

ans(16)  =  ans(8)/s4 
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RETURN 
END 
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APPENDIX  B.  HENON  EQUATION  ANALYSIS  PROGRAM 

Q      **************************************************************** 

PROGRAM  CHAOS 
C 

C  Written  by  Ray  Kamada  &  Korey  Jackson 
C 

DOUBLE  PRECISION  X(5000),  Y(5000),  DA,  A,  B, 
+   DAERR,  L(100),  Bl,  S 

INTEGER  NMAX, II, JMAX, I FLAG,  Jl,  NMAX1,  COUNT,  EPSIL(IOO) 

OPEN(2,  FILE=' LENGTH. DAT' ,  STATUS='NEW,  ERR=150) 

WRITE (2,*)'    A      B  DA  DAERR' 

OPEN(7,  FILE=' ENTROP.DAT' ,  STATUS= 'NEW* ) 

WRITE(7,*)'    A      B  S' 

C 

C  Establish  NMAX 
C 

NMAX=3600+1 

NMAX1=NMAX+100 

JMAX=NMAX/2 

PRINT* , ' JMAX= ' , JMAX 

IF  ( NMAX. GT. 5000)  GO  TO  100 
C  Initialize  all  variables 
IFLAG=0 

CALL  DIVISR(NMAX,  JMAX,  EPSIL,  COUNT) 
C 

C  Now  set  up  a  chaos  generating  routine. . . 
C 

A=1.0570D0 

B1=0.3D0 
11=0 
C30     IF ( II. GT. 600) GOTO  55 

J1=0 
C  A=A+0.0010D0 

B=B1 

11=11  +  1 
c  35        IF(J1.GT.30)GOTO  50 
c  J1=J1+1 

c  B=B+0.05D0 

print*, 'MAIN1   A,  B,',A,B 
IFLAG=0 

CALL  HENON(NMAXl,A,B,X,Y,IFLAG) 

print*, 'henon  complete' ,a,b, if lag 
C  Error  trap  for  x,  y  out  of  range  (diverging  solutions) 
c        IF(IFLAG.EQ.l)GOTO  35 

C  Comment  out  above  line  &  replace  w/below  line  when  B  is  fixed 
C       IF(IFLAG.EQ.l)GOTO  30 

if (if lag.eq. 1. )goto  200 
C 

C  Analyze  the  data 
C 

CALL  ANALYS(NMAX, JMAX,X, Y, A, B,  DA,  DAERR,  COUNT,  EPSIL) 

print* ,' analysis  complete',  a,  b,  DA 

44  WRITE (2, 45)  A,  B,  DA,  DAERR 

45  FORMAT(  F7.4,  3X,  F6.3,  5X,  F8 . 5 ,  5X,  F8.5) 

46  print*, 'WRITE  to  file  complete  in  analysis' 
CALL  ENTROP(NMAX,X,Y,S) 
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WRITE(7,47)  A,  B,  S 
47     FORMAT (F7. 4, 3X,  F6. 3 , 5X, F8 . 5 ) 
c  49     GOTO  35 
50     CONTINUE 
C  54     GOTO  30 
C 

55  CONTINUE 
C  CLOSE ( 2 ) 
C       CLOSE ( 7 ) 

GOTO  200 
C 
100    PRINT*, 'ERROR NMAX  Greater  than  10000 Array  size  too  small' 

GOTO  200 
146    PRINT*, 'ERROR. . .Can  not  open  henon.dat  file' 

GOTO  200 
148    PRINT*, 'ERROR. . .Can  not  open  analys.dat  file' 

GOTO  200 
150   PRINT*, 'ERROR. . .Can  not  open  length.dat  file' 
200    PRINT*, 'COMPLETED  MAIN2 ' 

END 

Q      **************************************************** 

SUBROUTINE  DIVISR(NMAX, JMAX,EPSIL,  COUNT) 

INTEGER  CHECK,  CHECK1,  COUNT,  NMAX,  JMAX,  EPSIL(IOO) 

CHECK=0 

CHECK1=0 

COUNT=0 

NMAX=NMAX-1 
1000   CHECK=CHECK+1 

IF ( CHECK. GT. JMAX)  GOTO  1050 

CHECK1=NMAX/CHECK 

CHECK1=CHECK1*CHECK 

IF(CHECKl.NE.NMAX)  GOTO  1000 

COUNT=COUNT+l 

EPSIL ( COUNT ) =CHECK 

PRINT*, 'EPSIL=' ,EPSIL(COUNT) ,  COUNT 

GOTO  1000 
1050   CONTINUE 

NMAX=NMAX+1 

RETURN 

END 

Q      ***************************************************************** 

SUBROUTINE  HENON(NMAXl, A, B,X, Y, IFLAG) 

DOUBLE  PRECISION  X(NMAXl),  Y(NMAXl),  A,  B,  CHECK 

INTEGER  NMAX1,  N,  IFLAG,  NMAX 
C      OPEN ( 4, FILE=' HENON.DAT' ,  STATUS= ' NEW ,    ERR=700) 
C 

C   Initialize  variables 
C 

PRINT*, 'PASSED  TO  HENON  OK' 

DO  300  N=1,NMAX1 
X(N)=0.0D0 
Y(N)=0.0D0 
300    CONTINUE 
C 

C  MAIN  COMPUTATION  ROUTINE 
C 

C   First  initiate  X(l),  if  desired. 
X(l)=-0.65d0 
Y(l)=  0.38d0 
C 
c      WRITE (4,*)       N  X  Y' 
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DO  400  N=l,  NMAX1 

X(N+1)=1.-A*X(N)*X(N)+Y(N) 
CHECK=X(N+1) 

IF ( CHECK. GE. 100. OR. CHECK. LE. -100. )  THEN 
IFLAG=1 

PRINT*, 'PROGRAM  OUT  OF  RANGE  IN  X' 
ENDIF 

Y(N+1)=B*X(N) 
CHECK=Y(N+1) 

IF ( CHECK. GE . 1000 . OR . CHECK. LE . -1000 . 0 )  THEN 
IFLAG=1 

PRINT*, 'PROGRAM  OUT  OF  RANGE  IN  Y' 
ENDIF 

IF(IFLAG.EQ.l)GOTO  760 
IF ( N . LT . 10 1 ) GOTO400 
C      WRITE(4,350)  N,  X(N) ,  Y(N) 
350    FORMAT(I10,5X,E15.5,5X,E15.5) 
400    CONTINUE 
C 

C   NOW  SHIFT  THE  VALUES  DOWN  IN  X(Y),  Y(I) 
C 

NMAX=NMAX1-100 
DO  600  N=l,  NMAX 
X(N)=X(N+100) 
Y(N)=Y(N+100) 
600     CONTINUE 
C 

C  Return  to  the  main  program 
C 

PRINT*, 'COMPLETED  HENON  OK  FOR',A,B 
GO  TO  770 
C 

700    PRINT*, 'Error  in  opening  HENON.DAT  file' 
C  740    CLOSE(2) 
C  750    CONTINUE 

GOTO  770 
760    PRINT*, 'X  value  out  of  range... A, B  TOSSED',  A,  B 

IFLAG=1 
770    CLOSE (1) 
RETURN 
END 

Q      It**************************************************************** 

SUBROUTINE  ANALYS (NMAX, JMAX,  X,Y,A,B,  DA,  DAERR,  COUNT, 
+  EPSIL) 

DOUBLE  PRECISION  SCRAP1,  SCRAP2 ,  L(100),  X(NMAX),  SUM (100), 
+   Y(NMAX),  EPS  (100),  DA,  A-,  B,  DAERR,  SCRAP3,  XY(200), 
+   ANS(16) 

INTEGER  NMAX,  I,  J,  CHECK,  CHE CK 1 , JMAX , COUNT ,  EPSIL (100), 
+   COUNT2 ,  nyval ( 100 ) ,  COUNT3 
C 

C    VARIABLE  LIST: 

C  NMAX=MAXIMUM  NUMBER  OF  DATA  POINTS 

C  COUNT=NUMBER  OF  EVEN  DATA  POINTS  THAT  WILL 

C  DIVIDE  EVENLY  INTO  NMAX 

C  L(J)=LENGTH  FOR  A  GIVEN  EPSILON 

C  EPSIL(J)=WIDTH,  OR  AN  INTEGER  DIVISOR  OF  NMAX 

C  EPS(J)=EPSIL(J)  IN  INTEGER  FORM 

C  DAERR=STANDARD  DEVIATION  OF  DA  (LOG  FORM) 

C 

OPEN  (3,  FILE= 'ANALYS. DAT ' ,  STATUS= ' NEW , 
+  ERR=9010) 
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c 

C      INITITALIZE  THE  VARIABLES 
C 

DO  4000  J=l,  100 
L(J)=0.0D0 
SUM(J)=0.0D0 
4000    CONTINUE 
C 

C     Now  compute  L(epsil) 
C 

DO  7500  J=l,  COUNT 
CHECK=NMAX-EPSIL(J) 
CHECK1=EPSIL(J) 
DO  7000  I=CHECK1,  CHECK,  CHECK1 
SCRAP1=X ( 1+1 ) -X ( I-CHECK1+1 ) 
SCRAP2=Y ( 1+1 ) -Y ( I-CHECK1+1 ) 
SCRAP1=SCRAP1*SCRAP1 
SCRAP2=SCRAP2  *SCRAP2 
SCRAP2=SCRAP2+SCRAP1 
SCRAP2=DSQRT ( SCRAP2 ) 
SUM ( J ) =SUM ( J ) +SCRAP2 
7000   CONTINUE 

L(J)=SUM(J) 
7  500   CONTINUE 
C 
C   QUIK  DATA  PRINT 

DO  7550  J=l,  COUNT 
L(J)=L(J)+1.0D0 
SCRAP1=L(J) 
L(J)=  DLOGIO(SCRAPI) 
SCRAP3=DBLE(EPSIL(J) ) 
EPS ( J ) =DLOG10 ( SCRAP3 ) 
WRITE(3,7545)  EPS(J),  L(J) 
7545   FORMAT(5X,F8.5,5X,F8.5) 
7550   CONTINUE 
C 

C  Now  do  a  linear  regression  to  find  DA 

C   First  put  into  a  format  to  use  the  canned  linear  regression 
C     subroutine. 

COUNT3=2*COUNT 
DO  7700  J=2,  COUNT3,  2 
I=J/2 

XY(J-1)=EPS(I) 
XY(J)=L(I) 
7700   CONTINUE 
C 

c  COUNT2  IS  THE  COUNT  WITH  THE  LAST  3  VALUES  TOSSED  FOR  Da 
c    computation 

COUNT2=COUNT-3 

CALL  LINREG  (COUNT2,  nyval,  XY,  ANS) 
print*, 'made  it  out  of  linreg' 
DA=ANS ( 2 ) 
DA=-DA 

DAERR=ANS(14) 
GOTO  9100 
9000   PRINT*, 'ERROR  OPENING  HENON.DAT  FILE' 
9010   PRINT*, 'ERROR  OPENING  ANALYS.DAT  FILE' 
9020   PRINT*, 'ERROR  OPENING  LENGTH.DAT  FILE' 
9100   print*, 'made  it  to  9100' 
RETURN 
END 
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Q      ***************************************************************** 

subroutine  LINREG  (k,  nyval,  xy,  ans) 
c 

c  This  subroutine  computes  the  slope  and  other  statistics  for  a 

c  linear  regression  with  several  y  values  for  each  x  value  or  with 

c  one  independent  variable, 
c 

c  argument   use     description 

c  k       input   Number  of  different  x  values 
c 

c  nyval    input   Array  of  number  of  y  values  for  each  x 

c  value.  Dimensioned  k 

c  xy      input   array  of  x  and  y  pairs  (xl,  yl,  x2,  y2,  etc.) 

c  of  dimension  (2,k) 
c 

c  ans     output  Array  of  results  computed,  dimensioned  16. 
c 

c  ans 

c  1  Number 

c  *2  Slope 

c  3  Mean  of  y 

c  4  Mean  of  x 

c  5  y-intercept 

c  6  Sum  of  y**2 

c  7  Sum  of  squares  mean 

c  8  Sum  of  squares  slope 

c  9  Residual 

c  10  standard  deviation  of  x 

c  11  standard  deviation  of  y 

c  12  standard  deviation  error 

c  13  standard  deviation  y-bar 

c  *14  standard  deviation  slope 

c  15  standard  deviation  y-intercept 

c  16  F-Ratio  for  slope 

c  **************************************************************** 

DOUBLE  PRECISION  nyval ( 1 ) ,  xy(l),  ans (16),  si,  s2,  sumx,  sumx2 
+  sumy,  sumy2,  sumxy,  enl,  en2,  s4,  endsl,  al3,  al4 

sumx  =  0.0D0 

sumy  =  0.0D0 

sumx2  =  0.0D0 

ans (6)  =  0.0D0 

sumxy  =  0.0D0 

n  =  0 

nx  =  1 

ny  =  2 

DO  20  j  =  l,k 

nyval (j)  =  1 

m  =  nyval ( j ) 

n  =  n  +  m 

DO  10  i  =  1,  m 

sumx  =  sumx  +  xy(nx) 

sumx2  =  sumx2  +  xy (nx) *xy (nx) 

sumy  =  sumy  +  xy(ny) 

ans (6)  =  ans (6)  +  xy (ny ) *xy (ny ) 

sumxy  =  sumxy  +  xy (nx) *xy (ny ) 
10     ny  =  ny  +  1 

nx  =  nx  +  nyval ( j )  +1 
20     ny  =  ny  +  1 

en  =  n 

ans(l)  =  en 

si  =  ans(l)*sumx2  -  sumx* sumx 
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s2  =  ans(l)*sumxy  -  sumx*sumy 

enl  =  en  -  1.0 

en2  =  enl  -  1.0 

ans(2)  =  s2/sl 

ans(3)  =  sumy/en 

ans(4)  =  sumx/en 

ans(5)  =  ans(3)  -  ans (2 ) *ans(4) 

ans ( 7 )  =  ans ( 3 ) *  sumy 

ans(8)  =  ans(2)*s2/en 

ans (9)  =  ans (6)  -  ans(7)  -  ans (8) 

s4  =  ans (9) /(en  -  2.0) 

endsl  =  en/sl 

ans (10)  =  sqrt(1.0/(enl*endsl) ) 

ans (11)  =  sqrt((ans(6)  -  ans(7))/enl) 

ans (12)  =  sqrt(s4) 

al3  =  s4/en 

ans (13)  =  sqrt(al3) 

al4  =  s4*endsl 

ans (14)  =  sqrt(al4) 

ans (15)  =  sqrt(al3  +  al4*ans (4) *ans(4)  ) 

ans(16)  =  ans(8)/s4 

RETURN 

END 

Q      **************************************************************** 

C   23  Apr  92  revision  (goes  with  0<a<1.25  plot) 
C 

SUBROUTINE  ENTROP ( NMAX , X , Y, S) 

DOUBLE  PRECISION  X(NMAX),  Y(NMAX),  S,  XMAX,  XMIN,  YMAX,  YMIN, 
+   PROB(200,200),  INVNMA,  P,  XRANGE,  YRANGE 
INTEGER  IX,  IY,  I,  J,  NMAX 
C  First  find  maximum  and  minimum  values 
XMAX=0.0D0 
YMAX=0.0D0 
XMIN=0.0D0 
YMIN=0.0D0 

print*, ' finding  max  and  minimum  values' 
DO  2000  1=1,  NMAX 

IF(X(I) .GT.XMAX)  XMAX=X(I) 
IF(X(I) .LT.XMIN)  XMIN=X(I) 
IF(Y(I) .GT.YMAX)  YMAX=Y ( I ) 
IF(Y(I) .LT.YMIN)  YMIN=Y(I) 
2000   CONTINUE 

C  Now  initialize  the  probability  array 
print*, ' initializing  P  array' 
DO  2105  1=1,  200 

DO  2100  J=l,  200 
PROB(I,J)=0.0D0 
2100      CONTINUE 
2105   CONTINUE 

C  Count  the  number  of  points  in  each  unit  area 
C  First  set  up  the  scheme 
XRANGE=XMAX-XMIN 
XRANGE=2 . 00D2 /XRANGE 
YRANGE=YMAX-YMIN 
YRANGE=2 . 00D2 /YRANGE 
DO  2110  1=1,  NMAX 

IX=INT ( ( X ( I ) -XMIN ) *XRANGE ) +1 
IY=INT( (Y( I) -YMIN)* YRANGE )+l 
IF(IX.GT.200)  IX=200 
IF(IY.GT.200)  IY=200 
PROB ( IX, IY ) =PROB ( IX, IY ) +1 . 0D0 
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2110   CONTINUE 
C  Now  normalize 

INVNMA=1 . 0D00/DBLE ( NMAX ) 
DO  2120  IX=1,  200 

DO  2115  IY=1,  200 

PROB(IX,  IY)=PROB(IX,IY)*INVNMA 
2115      CONTINUE 
2120   CONTINUE 

C  Finally,  compute  the  entropy  S 
S=0.0D0 

print*, ' computing  S' 
DO  2130  IX=1,  200 

DO  2125  IY=1,  200 

IF(P.ne.0d0)print*, 'P=' ,P 
P=PROB(IX,IY) 
IF(P.EQ.ODO)  GOTO  2125 
S=S-P*DLOG(P) 
2125      CONTINUE 
2130   CONTINUE 

C  Now  return  to  the  main  program  to  write  S  for  this  value  of  A 
RETURN 
END 


132 


APPENDIX  C.  PROGRAM  MCNIDER 

************************************************************************ 

PROGRAM  MCNIDE 
c 

c  This  program 

c  1)  computes  Monte  Carlo  particle  velocity  and  position  as  a  function 

c  of  time  step,  t+A*t,  using  the  McNider  algorithm,  Kamada  vertical 

c  velocity  skewness  variations,  and  Kamada  mesoscale  windflow  and 

c  turbulence  simulation  model, 
c 

c  2)  computes  the  chaos  metrics,  Da  (self-affine  fractal  dimension) 

c  for  velocity,  S  (information  entropy),  and  lambda  (the  Lyapunov 

c  exponent),  as  well  as  Ri,  the  atmospheric  Richardson  number,  BVF 

c  ( Brunt -Vaisala  frequency)  1  (buoyancy  length  scale),  au,v,w,  the 

c  velocity  variances,  TKE,  the  turbulence  kinetic  energy,  Db 

c  (self-similar  fractal  dimension  for  position). 


& 

St 

& 
OPEN 
OPEN 
OPEN 
OPEN 


el,  p(1000),  entropy ( 101) ,   up,  dn,  lyap,  pu,  pv,  pw, 
sum(45,  100),  DavsL(8,100) ,  pk,  ans(16),  xy(90) 
'ALENCE  (vert,  r  ) 

k/1,  2,  3,  4,  5,  6,  8,  9,  10,  12,  15,  16,  18,  20,  24,  25, 
30,  36,  40,  45,  48,  50,  60,  72,  75,  80,  90,  100, 
.20,  144,  150,  180,  200,  225,  240,  300,  360,  400, 
50,  600,  720,  900,  1200,  1800,  3600/ 

(1,  FILE=' McNideis.dat' ,  RECL=255,  STATUS= ■ new' ,  ERR=700) 
(2,  FILE=' McNidell.dat' ,  RECL=255,  STATUS= ' new' ,  ERR=800) 
(3,  FILE='DavsL.daf ,  RECL=255,  STATUS= ' new ' ,  ERR=900) 
OPEN  (4,  FILE='SvsL.daf ,  RECL=255,  STATUS= ' new' ,  ERR=1000) 
OPEN  (7,  FILE=*lyap.daf ,  RECL=255,  STATUS= ' new' ,  ERR=1100) 
OPEN  (8,  FILE=' bifurct.dat' ,  RECL=255,  STATUS= ' new ' ,  ERR=1200) 
OPEN  (9,  FILE=' initial. daf  ,  RECL=255,  STATUS= ' old' ,  ERR=1300) 
c 
c  ********** 

c   INITIALIZE 
c   ********** 

READ  (9,30) 
30     FORMAT  (///) 

READ  (9,*)  iis,  ill,  iDavsL,  isvsL,  ilyap,  ibifurct,  noskew, 
&    iMcNid,  iKamada,  dt,  ivert,  iphase,  iveloc,  windspd2,  zO, 
&   Start,  m2 ,  iw,  ix,  upf actor,  dnf actor,  dwf actor,  zinitial, 
St    icount,  iposi 

WRITE ( 6, * ) ' iis,  ill,  iDavsL,  isvsL,  ilyap,  ibifurct,  noskew,' 

WRITE (6,*) 'iMcNid,  iKamada, dt, ivert , iphase, iveloc,  windspd2,z0, ' 

WRITE(6, *) ' Start,  m2,  iw,  ix,upfactor,  dnf actor, dwf actor, zinitial ' 

WRITE ( 6, *)' icount,  iposi' 

WRITE (6,*)  iis,  ill,  iDavsL,  isvsL,  ilyap,  ibifurct,  noskew, 
St    iMcNid,  iKamada,  dt,  ivert,  iphase,  iveloc,  windspd2 ,  zO, 
St   Start,  m2 ,  iw,  ix,  upf  actor,  dnf  actor,  dwf  actor,  zinitial, 
St   icount,  iposi 
c   Comments  on  input  parameters: 
c         iis  "     =  1  switch  gives  position  and/or  velocity  file 
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iil       =  1  switch  gives  Ln  e  vs  Ln(L(e)  file 

iDavsL    =  1  switch  gives  Da  vs  1/L  file 

isvsL     =  1  switch  gives  entropy,  S,  vs  1/L  file 

ilyap    =  1  switch  gives  Lyapunov  exponent  vs  1/L  file 

ibifurct  =  1  switch  gives  bifurcation  disagram  file 

noskew   =  1  switch  turns  off  vertical  velocity  skewness 

computation 
iMcNid   =  1  switch  turns  on  McNider  vertical  velocity  skewness 

scheme 
iKamada  =  1  switch  turns  on  Kamada  vertical  velocity  skewness 

scheme 
dt       =  time  step  size  in  seconds 
ivert    =  1  switch  gives  vertical  height  file 
iphase   =  1  switch  gives  phase  space  vector  file 
iveloc   =  1  switch  gives  velocity  vector  file 
windspd2  =  input  proscribed  wind  speed  at  2  meters 
zi       =  input  proscribed  fully  turbulent  boundary  layer 

height 
zO       =  input  proscribed  surface  roughness  height  ("1/7)  the 
average  height  of  surface  roughness  elements  if 

element 
density  is  between  10  and  50%  (like  most  vegetation) 

Do  some  idiot  proofing 

IF  (  ivert  .eq.  1)  THEN 

iphase  =  0 

iveloc  =  0 

iposi   =  0 

END  IF 
IF  (  iphase 

ivert  =  0 

iveloc  =  0 

iposi   =  0 

END  IF 
IF  (  iposi  .eq 

ivert  =  0 

iphase  =  0 

iveloc  =  0 

ENDIF 
IF  (  noskew 

iMcNid  =  0 

iKamada=  0 

ENDIF 
IF  (  iMcNid 

noskew  =  0 

iKamada=  0 

ENDIF 
IF  (  iKamada. eq. 

noskew  =  0 

iMcNid  =  0 

ENDIF 
c   Initialize  variables 
ml  =  1 
n     =  3600 
en   =  n 

WRITE  (4,*)'      S 
WRITE  (7,*)'     Lyap 
WRITE  (8,*)  '    1/L 
DO  500  1  =  ml,  m2 
c   Let  Linv,  the  inverse  Obukhov  length,  1/L  =  -kgw'O' / (0*ustar**3 ) ,  be 


eq.  1)  THEN 


1)  THEN 


eq.  1)  THEN 


eq.  1)  THEN 


1)  THEN 


1/L' 

1/L' 
bifurcation  pt ' 
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c   varied  in  increments  along  the  interval  from  -1  to  1. 
em2  =  m2 

Linv  =  start  +  0.3d0* ( 1-1) /em2 
zi  =  2000  -  6000* (Linv  +  .2) 
r(0)   =  0. 
rl(0)  =  0. 
r2(0)  =  0. 
vert(0)  =  zinitial 
x  =  0.001 
y  =  0.002 
z  =  zinitial 
udp  =  0.001 
vdp  =  0.002 
wdp  =  0.003 

Q  ******** 

c  GET  iteration  series  for  a  given  1/L 
c   ******** 

IF  (  iis  .ne.  1  )  GOTO  150 
pu  =  0.0 
pv  =  0.001 
c        pw  =  -.01*alogl0(l) 
pw  =  -0.001 
time  =0.0 
rmax  =  0.0 
r2max  =  0.0 
rmin  =  0.0 
r2min  =  0.0 
sigU2sum  =  0. 
sigV2sum  =  0. 
sigW2sum  =  0. 
zsum  =  0. 
BVFsum  =  0. 
BLSsum  =  0. 
sumdw  =  0 . 
tkebar  =  0. 
wbuoy  =  0. 
sumw  =  0 . 
sumwbuoy  =  0 . 
d        WRITE (6,*) ' rmax2  =',rmax2 
rmax2  =  0 . 
rmin2  =  0. 
iflag  =  0 
DO  100  i  =  1,  n  +  100 

CALL  McNid  (x,  y,  z,  udp,  vdp,  wdp,  dt,  noskew,  iMcNid, 
£r      iKamada,  Linv,  windspd2 ,  deltaz,  zi,  zO,  wbuoy,  sumw,  Ri, 
&       sumwbuoy,  BVF,  BLS,  sigmau2,  sigmav2,  sigmaw2,  tke,  icount, 
&      pu,  pv,  pw,  i,  time,  wk,  ustar,  iw,  ix,  upfactor,  iflag, 
&       sigmau,  sigmav,  sigmaw,  Tlw,  sumdw,  dnfactor,  dwfactor  ) 
BVFsum  =  BVFsum  +  BVF 
BLSsum  =  BLSsum  +  BLS 
sigU2sum  =  sigU2sum  +  sigmau**2 
sigV2sum  =  sigV2sum  +  sigmav**2 
sigW2sum  =  sigW2sum  +  sigmaw**2 
zsum  =  zsum  +  z 

bltime  =  0. 3*zi/ ( sigmaw*3700. ) 
eddytime  =  0.01*Tlw 
c  IF  (Linv  .ge.  0.000)  dt  =  MIN  (eddytime,  bltime) 

c  IF  (Linv  .ge.  0.000)  dt  =  MIN  (dt,  0.25) 

c  IF  (Linv  .gt.  0.000)  dt  =  eddytime 

c  IF  (Linv  .It.  -0.00)  dt  =  MIN  ( 1. 2*zi/ ( 3700. 0*wk) ,  eddytime) 

dt  =  0.166666667 
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time  =  time  +  dt 
cd  WRITE  (6,*)  'time', time 

c   get  position  vector,  rl 

rl(i)  =  sqrt(  x*x  +  y*y 
c   get  velocity  vector,  r2 

r2(i)  =  sqrt(  udp*dup  + 
c   get  vertical  position  series 

IF  (  ivert  .eq.  1)  r(i) 
c   get  position  series 

IF  (  iposi  .eq.  1)  r(i) 
c  get  phase  space  position 

IF  (iphase.eq.  1)  r(i)  = 
c   get  velocity  series 

IF  (iveloc  .eq.  1)  r(i) 
c   get  min/max  range  of  the  values 

IF  (  i  .gt.  100  .and.  r(i) 

IF  (  i  .gt.  100  .and.  r(i) 

IF  (  i  .gt.  100  .and.  r(i) 

IF  (  i  .gt.  100  .and.  r(i) 
c  WRITE  (1,50)  i, 


+  (z-zinitial)**2 
vdp*vdp  +  wdp*wdp 
=  z  -  zinitial 
-  rl(i) 

sqrt(  rl(i)*rl(i) 
=  r2(i) 

) 


+  r2(i)*r2(i)  ) 


.gt. 
.It. 
.gt . 
•  It. 


c50 
100 
150 


14, 


rl(i) 
3(3x, 


r2(i) 
fl2.5) 


rmax 
rmin  ) 
rmax2 ) 
rmin2 ) 
r(i) 


rmax  = 
rmin  = 
rmax 2= 
rmin2= 


r(i) 
r(i) 
r(i) 
r(i) 


) 


Linv 


FORMAT  (lx, 

CONTINUE 

CONTINUE 
WRITE  (6, *)' iteration  series  created  for  1/L  = 
IF  (  ibifurct  „ne.  1  )  GOTO  240 

Q  *************** 

c   GET  BIFURCATION  diagram,  using  the  last  1600  points.   For  each  1/L 
value 

points  in  "bi" 


put  all  the  unique 
*************** 


m  =  1 
bi(l,m)  = 
bi(2,m)  = 


Check 


Check 


Linv 
r(1999) 
slop  =  0.01* (rmax  -  rmin) 
IF  (iveloc  .eq.  1)  slop  = 
eq.  1)  slop  = 
eq.  1)  slop  = 


IF  (ivert 
IF  ( iphase 
every  r(i) 
DO  220  i  = 
each  r(i+l) 
DO  180  j 


0.005*11. 

0.005*zi 

0.002*3000. 


2000,  n+99 

against  all 
=  1,  m 


unique  points  in  bi 


c 
c 
180 


220 
c 

230 
240 


up 

dn 

IF 

Since  r(i+l) 

unique  point, 


=  bi(2,j)  +  slop 
=  bi(2,j)  -  slop 
(  r(i+l)  .gt.  dn  .and. 
lies  beyond  margins  of 
so  add  it  to  the  list 


r(i+l)  .It.  up  )  GOTO  220 
all  "bi",  we  assume  it  a  new 
of  values  in  "bi" . 


CONTINUE 
bi(2,m+l)  =  r(i+l) 


bi(l,m+l)  =  Linv 

m  =  m+1 

CONTINUE 
IF  (  iveloc  .ne.  1  )  bi(2,m)  =  r(1999) 
WRITE  (8,230)  ((  bi(i,j),  i  =  1,2),  j  =  l,m) 
FORMAT  (  lx,  f9.5,  5x,  f 9 . 5  ) 
CONTINUE 

WRITE  (6, *) 'bifurcation  diagram  points  found 
IF  (  ilyap  .ne.  1  )  GOTO  260 


for  1/L  = 


Linv 


************ 


GET  LYAPUNOV  exponent  for  this  series, 


************ 
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c  n 

c      lyap  =  lim  (1/n)  sum  log  |f'(x    )|   ,  where  f'(x    )  =^i  ( 1  -2x   ) 

c  n>»       i=i    e     i+1  i+1  i+1 

c 

c    where  we  start  from  y(100)  to  avoid  dependence  on  initial  condition, 

c     xO. 

c 

lyap  =  0.0 

DO  250  i  =  1,  n 
c  lyap  =  lyap  +  dlog(  abs(mu*(1.0  -  2.*y(i+100)  )  )) 

lyap  =  0.0000000000 
250  CONTINUE 

c    Normalize  the  lyapunov  exponent  by  N  -  100,  the  number  of  points 
c    scanned. 

lyap  =  lyap/en 

WRITE  (7,255)  lyap,  Linv 
255  FORMAT  (lx,  2(f9.5,2x)  ) 
2  60  CONTINUE 

IF  (  isvsL  .ne.  1  )  GOTO  290 
c 

C  ***********  ff 

c  GET  ENTROPY  value,  S  =  -  sum  p  log  p   ,  for  this  iteration  series  &  1/L 
c  ***********  i=l   i   e  i 

c 

c   Zero  out  the  number  of  points  in  each  bin 
DO  265  ii  =  1,  100 
ipts(ii)  =  0 
265  CONTINUE 

c     Subdivide  the  interval  into  100  sub-intervals, 
c     ii,  and  check  all  pts(i),  from  101  -  3600. 
entropy (1)  =  0.0 
DO  270  i  =  101,  n  +  100 

ii  =  99* (  r(i)  -  rmin2 ) / ( rmax2  -  rmin2 )  +  1.01 
c     Make  sure  we  don't  create  extra  sub-intervals 

IF  (ii  .gt.  100)  ii  -  100 
c     Increment  the  ipts  counter  for  sub-interval  ii  for  every  r(i)  found 
c     in  ii. 

ipts(ii)  =  ipts(ii)  +  1 
270  CONTINUE 

DO  285  ii  =  1,  100 
c     If  no  points  are  in  sub-interval,  ii,  leave  the  entropy  unchanged. 
IF  (  ipts(ii)  .eq.  0)   THEN 
si  =  0.0 
GOTO  280 
ENDIF 
c     Otherwise  compute  probability  of  a  point  being  in  ii  as  the  #  of 
c     points  in  ii  divided  by  N,  the  total  number  of  points  scanned, 
pts  =  ipts(ii) 
prob  =  pts/en 
si  =  prob*alog(prob) 
c     Since  the  probability  is  less  than  or  equal  to  1,  the  log  is 
c     negative,  so  si  is  less  than  or  equal  to  zero,  so  the  entropy  will 
c     grow  more  positive,  if  the  r(i)  are  scattered  over  more  than  one 
c     sub-interval. 

280         entropy (1)  =  entropy (1)  -  si 
285  CONTINUE 

WRITE  (4,255)  entropy (1),  Linv 

IF  (  ill.  eq.  1)  WRITE  (2,  *)'    Ln(e)    Ln  [L(e)]' 
WRITE  (6, *) 'entropy  value  found  for  1/L  =  ',  Linv 
290       CONTINUE 

IF  (  iDavsL  .ne.  1  )  GOTO  500 
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Q        ******* 

c  GET  Da.  First  get  total  length  as  sum  of  y(i)  changes  for  varying 
c  *******  spacings  of  i. 
nn  =  45 
DO  400  j  =  1,  nn 

sum(j,l)  =  1.0000 

DO  300  i  =  k(j),  3600,  k(j) 

sum(j,l)  =  sum(j,l)  +  abs (  r(i+100)  -  r(i-k(j)+100  )) 
300  CONTINUE 

c    Get  log  of  the  total  length 

lglngth(j,l)  =  dlogl0(  sum(j,l)  ) 
c    Get  log  of  the  spacing 
pk  =  k(j) 

lgk(j)   -  dlogl0(  pk  ) 
c    Put  these  log  values  into  an  array  and  send  to  the  linear  regression 
c    routine. 

xy(2*j-l)  =  lgk(j) 
xy(2*j)  =  lglngth(j,l) 

IF  (ill  .eq.  1)    WRITE  (2,350)  lgk(j),  lglngth(j,l) 
350  FORMAT  (lx,  f8.5,  3x,  f8.5) 

400  CONTINUE 

c    WRITE  (6,*)  'program  passed  statement  400' 
c 

c    DO  LINEAR  REGRESSION  to  get  Da  from  loglO(k)  vs.  loglO  (L(e)]  values 
c 

c    Use  only  spacings  up  to  1200  points  for  the  linear  regression.   So 
c    skip  spacings  1800  and  3600. 
nn  =  nn  -  2 

CALL  STLNRG(nn,  nyval,  xy,  ans) 
c 

c   -ans (2)  contains  Da 

c   ans (14)  contains  the  standard  deviation  of  Da 
c 

c   Put  1/L,  Da,  aDa,  BVF,  BLS,  aw1,  tke  in  a  7  by  1  array 
c 

ptsinv  =  l./(n+100.) 
BVFbar  =  BVFsum*ptsinv 
BLSbar  =  BLSsum*ptsinv 
sigw2bar  =  sigw2sum*ptsinv 

tkebar  =  . 5*(sigu2sum  +  sigv2sum  +  s igw2 sum )* ptsinv 
zbar  =  zsum*ptsinv 
DavsL(l,l)  =  -ans (2) 
DavsL(2,l)  =  -ans (14) 
DavsL(3,l)  =  Linv 
DavsL(4,l)  =  BVFbar 
DavsL(5,l)  =  BLSbar 
DavsL(6,l)  =  sigw2bar 
DavsL(7,l)  =  tkebar 
DavsL(8,l)  =  zbar 

WRITE  (6,*) 'Da  value  found  for  1/L  =  ',  Linv 
500       CONTINUE 

WRITE  (3,*)'     Da        cDa         Linv       BVF       BLS' 
&  , '        aw2    tke    zavg ' 

WRITE  (6,*)'     Da         oDa  Linv        BVF        BLS ' 

&  , '        ow2    tke    zavg ' 

WRITE  (3,550)  ( (DavsL( i, 1 ) ,  i  =  1,8),  1  =  ml,m2) 
WRITE  (6,550)  ( (DavsL( i, 1) ,  i  =  1,8),  1  =  ml,m2) 
550    FORMAT  (lx,  5(fl0.4,lx),  2f7.2,  f7.0) 

STOP 
700    WRITE  (6,*) 'error  opening  McNideis.dat' 
STOP 
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800    WRITE  (6,*)'error  opening  McNidell.dat' 

STOP 
900   WRITE  (6,*)' error  opening  DavsL.dat ' 

STOP 
1000  WRITE  (6,*)'error  opening  SvsL.daf 

STOP 
1100  WRITE  (6,*)' error  opening  Lyap.daf 

STOP 
1200  WRITE  (6,*)' error  opening  bifurct.dat' 

STOP 
1300  WRITE  (6,*)'error  opening  initial.dat* 

STOP 

END 
*********************************************************************** 

SUBROUTINE  STLNRG  (k,  nyval,  xy,  ans) 
c 

c  This  subroutine  computes  the  slope  and  other  statistics  for  a 

c  linear  regression  with  several  y  values  for  each  x  value  or  with 

c  one  independent  variable, 
c 

c  argument   use     description 

c  k       input   Number  of  different  x  values 
c 

c  nyval    input   Array  of  number  of  y  values  for  each  x 

c  value.  Dimensioned  k 

c  xy      input    array  of  x  and  y  pairs  (xl,  yl,  x2,  y2,  etc.) 

c  of  dimension  (2,k) 
c 

c  ans     output  Array  of  results  computed,  dimensioned  16. 
c 

c  ans 

c  1  Number 

c  *2  Slope 

c  3  Mean  of  y 

c  4  Mean  of  x 

c  5  y- intercept 

c  6  Sum  of  y**2 

c  7  Sum  of  squares  mean 

c  8  Sum  of  squares  slope 

c  9  Residual 

c  10  standard  deviation  of  x 

c  11  standard  deviation  of  y 

c  12  standard  deviation  error 

c  13  standard  deviation  y-bar 

c  *14  standard  deviation  slope 

c  15  standard  deviation  y-intercept 

c  16  F-Ratio  for  slope 
c 
********************************************************************** 

DIMENSION  nyval (1) 
DOUBLE  PRECISION  xy ( 1 ) ,  ans (16) 
c     WRITE  (6,*)'  program  entered  linear  regression  routine' 
sumx  =  0.0 
sumy  =  0.0 
sumx2  =  0.0 
ans (6)  =0.0 
sumxy  =  0.0 
n  =  0 
nx  =  1 
ny  =  2 
DO  20  j  =  l,k 
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nyval(j)  =  1 

m  =  nyval ( j ) 

n  =  n  +  m 

DO  10  i  =  1,  m 

sumx  =  sumx  +  xy(nx) 

sumx2  =  sumx2  +  xy (nx) *xy (nx) 

sumy  =  sumy  +  xy(ny) 

ans(6)  =  ans(6)  +  xy (ny) *xy (ny) 

sumxy  =  sumxy  +  xy (nx) *xy (ny) 
10    ny  =  ny  +  1 

nx  =  nx  +  nyval ( j )  +1 
20    ny  =  ny  +  1 

en  =  n 

ans(l)  =  en 

si  =  ans(l)*sumx2  -  sumx*sumx 

s2  =  ans(l)*sumxy  -  sumx*sumy 

enl  =  en  -  1.0 

en2  =  enl  -  1.0 

ans(2)  =  s2/sl 

ans(3)  =  sumy/en 

ans (4)  =  sumx/en 

ans (5)  =  ans(3)  -  ans (2 ) *ans(4) 

ans (7)  =  ans(3)*sumy 

ans(8)  =  ans(2)*s2/en 

ans (9)  =  ans (6)  -  ans(7)  -  ans (8) 

s4  =  ans (9) /(en  -  2.0) 

endsl  =  en/sl 

ans (10)  =  sqrt(1.0/(enl*endsl) ) 

ans (11)  =  sqrt ( ( ans ( 6 )  -  ans(7))/enl) 

ans (12)  =  sqrt(s4) 

al3  =  s4/en 

ans (13)  =  sqrt(al3) 

al4  =  s4*endsl 

ans (14)  =  sqrt(al4) 

ans (15)  =  sqrt(al3  +  al4*ans(4) *ans (4)  ) 

ans(16)  =  ans(8)/s4 

RETURN 

END 
*********************************************************** 

SUBROUTINE  McNid  (x,  y,  z,  udp,  vdp,  wdp,  dt,  noskew,  iMcNid, 
&  IKamada,  Linv,  windspd2 ,  deltaz,  zi,  zO,  wbuoy,  sumw,  Ri, 
&  sumwbuoy,  BVF,  BLS,  sigmau2,  sigmav2,  sigmaw2,  tke,  icount, 
&  pu,  pv,  pw,  i,  time,  wk,  ustar,  iw,  ix,  upfactor,  iflag, 
&  sigmaudp,  sigmavdp,  sigmawdp,  Tlw,  sumdw,  dnfactor,  dwfactor) 

c  this  subroutine  computes  particle  velocity  and  position  for  different 

c  values  of  Ri#,  etc. 

c 

*********************************************************************** 

REAL  Km,  L,  lmdamu,  lmdamv,  lmdamw 

DOUBLE  PRECISION  Linv,  ruu,  rvv,  rww,  pu,  pv,  pw 
c  Get  mesoscale  and  boundary  layer  flow  parameters 
c  Initialize  variables 

zi  =  2000  -  6000* (Linv  +  .2) 

zovrzi  =  z/zi 

zovrL   =  z*Linv 

ziovrL  =  zi*Linv 

abszovrL  =  abs(zovrL) 

L  =  l./Linv 

CALL  MESOFLOW(  z,  zi,  zO,  Linv,  windspd2 ,  ustar,  Km,  Ri, 
&      u,  v,  w,  Theta2,  dudz,  dvdz,  wk,  D,  R,  wptpZ,  BVF,  iw, 
&     BLS,  sigmaU2,  sigmaV2,  sigmaW2,  tke,  Wz,  i,  wk,  ustar, 
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&     bdthetdz,Thstar,  Deltheta) 
d      IF( i. It . iw) WRITE (6, * ) ' zovrzi, zovrL, ziovrL' , zovrzi, zovrL, ziovrL 
c  Get  A,  lambdas,  ou,v,wdp,  and  buoyancy  length  scale,  BLS 
IF  (  zovrL  .It.  0.0  )  THEN 

sigmaudp  =  ustar*(12.0  +  0. 5*zi/abs (L) ) **0. 3333333 
d  IF(i.lt.iw)WRITE  ( 6, * ) '  cm  ' %=u* ( 12  +  . 5zi/ | L | ) **l/3 ' , 
d    &  sigmaudp, ustar, zi,L 

IF  (  abszovrL  .le.  1.  ) 
&     A  =  0.31*(1.  -  3.*zovrL)**(-0.333333)*(l.  -  15. *zovrL) **0.25* 
&  (0.55  +  0.38*zovrL  ) 

IF  (  0.1*abs( ziovrL)  .gt.  abszovrL  .and.  abszovrL  .gt.  1.) 
&     A  =  0.05*(1.  -  3.*zovrL)**(-0.333333)*(l.  -  15. *zovrL) **0. 25 
A  =  MAX  (  A,  0.07) 
A  =  MIN  (  A,  0.18) 
lmdamu  =  1.5*zi 

IF  (  z  .le.  0.1*zi  .and.  abs(zovrl)  .gt.  1.)  lmdamw  =  5.9*z 
IF  (  z  .le.  abs(L)  )  lmdamw  =  z/(  0.55  +  0.38*zovrL  ) 
IF  (  0.1*zi  .It.  z  .and.  z  .It.  zi  ) 
&      lmdamw  ■  1.8*zi*(  1.  -  exp(-4.0*zovrzi) 
&  -  0.0003*exp(  8.0*zovrzi  )  ) 

sigmawdp  =  Km/ (A*lmdamw) 
d        IF(i.lt.iw) 

d    &      WRITE (6,*) 'at  z,ow' ' =Km/( A* lmdamw) ' ,z, sigmawdp, Km, A, lmdamw 
ELSEIF  (  zovrL  .ge.  0.)  THEN 
sigmaudp  =  2.3*ustar 
lmdamu  =  0. 7* sqrt (zovrzi) *zi 
es  =  70.0 

IF  (z. It. 205.0)  es  =  0.35*z 
lmdamw  =  MAX(  z,  2.9*es) 
Ric  =  0.25 

Rifactor  =  ((Ric  -  1. 25*Ri) /Ric) **0. 58 
d         IF(i.lt.iw)WRITE(6,*) 'Rif=( (Ric-Ri) /Ric) ** . 58 ' , Rifactor, Ric, Ri 
shear  =  sqrt(dudz**2  +  dvdz**2) 
sigmawdp  =  1.2*es*Rif actor*shear 
d        IF(i.lt.iw)WRITE(6,*) 'ow' '=1.2esRi, shear ' , 
d    &  sigmawdp,  es,  Rifactor,  shear 

BLS  =  sigmawdp /BVF 
END  IF 
sigmavdp  =  sigmaudp 
IF  (ikamada  .eq.  1)  THEN 
c   Use  values  derived  from  Sorb j an  and  Kamada  from  mesoscale  subroutine 
d      IF(i.lt.iw)WRITE(6,*) ' ou2 ,ov2 , ow2 ' , sigmaU2, sigmaV2 , sigmaW2 
sigmaudp  =  sqrt (sigmaU2 ) 
sigmavdp  =  sqrt ( sigmaV2 ) 
sigmawdp  =  sqrt ( sigmaW2 ) 
d       IF( i. It . iw)WRITE(6, *)' iKmda  ou,v,w' ' ', sigmaudp, sigmavdp, sigmawdp 
c   To  get  vertical  integral  length,  see  pps.  244-247  Lectures  on  Air 
c   Pollution 

aln  =  0.41*z 

lmdamw  =  0.27*sigmawdp/BVF 
lmdamw  =  l./(l./aln  +  1. /lmdamw) 
d      IF(i. It. iw)WRITE(6,*) ' iKmda  aln, lmdamw' , aln, lmdamw 
ENDIF 

IF  (noskew  .eq.  1)  tke  =  0. 5* (sigmaU2  +  sigmaV2  +  sigmaW2) 
lmdamv  =  lmdamu 
d     IF ( i . It . iw ) WRITE ( 6 , * ) ' lmdamu , lmdamv , lmdamw ' , lmdamu , lmdamv , lmdamw 
c  Get  mean  velocity 

V  =  sqrt(  u*u  +  v*v  +  w*w  ) 
c  Get  betas 

d      IF(i.lt.iw)  WRITE  ( 6, *)' u,v,w, sigmaudp' , u,v,w, sigmaudp 
d      IF(i.lt.iw)  WRITE  ( 6, *)' sigmavdp, sigmawdp ', sigmavdp, sigmawdp 
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betau  =  0. 6*V/sigmaudp 

betav  =  0. 6*V/sigmavdp 

betaw  =  0. 6*V/sigmawdp 
d     IF(i.lt.iw)WRITE(6,*)  'at  z,/3w=.  6V/ow'  ' ',  z,  betaw,  V,  sigmawdp 
c  Get  Lagrangian  time  scales 
d     IF( i. It . iw)WRITE(6, *) 'betau, lmdamu, V , betau, lmdamu,V 

Tlu  =  0. 2 * betau* lmdamu /V 

Tlv  =  0.2*betav*lmdamv/V 

Tlw  =  0.2*betaw*lmdamw/V 
d      IF(i.lt.iw)WRITE(6,*)  'Tlw=. 2  *J3w*  lmdamw/ V  ,  z, Tlw, betaw,  lmdamw,V 

IF  (  ikamada  .eq.  1)  THEN 

IF  (  Linv  .ge.  0.  )  Tlw  =  lmdamw/ sigmawdp 
IF  (  Linv  .It.  0.  )  THEN 
Tlw  =  0.3*zi/wk 
AspectR  =  2.0  -  40.0*Linv 
Tlu  =  AspectR*Tlw 
Tlv  =  Tlu 
END  IF 
d      IF(i. It. iw)WRITE(6, *) ' Tlw= lmdamw/ ow' ' ' ,Tlw, lmdamw, sigmawdp 
END  IF 

IF  (  Tlu  .It.  0.02*dt  )  Tlu  =  0.02*dt 

IF  (  Tlv  .It.  0.02*dt  )  Tlv  =  0.02*dt 

IF  (  Tlw  .It.  0.02*dt  )  Tlw  =  0.02*dt 
c  Get  Autocorrelation  function 
d      IF(i.lt.iw)WRITE(6,*) ' dt ,TLu,TLv,TLw' ,  dt,  Tlu,Tlv,TLw 

Ru  =  exp(  -dt/Tlu) 

Rv  =  exp(  -dt/Tlv) 

Rw  =  exp(  -dt/Tlw) 
d     IF(i. It. iw) WRITE  (6,*)'at  z, Rw=exp( -dt/Tlw) ', z,Rw,dt, Tlw 
d     IF(i.lt. iw)WRITE  ( 6, *) 'Ru,Rv,Rw, ikamada' ,Ru,Rv,Rw, ikamada 
c  Get  standard  deviations  of  velocity 

sigmautp  =  sigmaudp*sqrt (  (1.  -  Ru**2)  ) 

sigmavtp  =  sigmavdp*sqrt (  (1.  -  Rv**2)  ) 

sigmawtp  =  sigmawdp*sqrt (  (1.  -  Rw**2)  ) 
d      IF(i.lt.iw) 

d    &   WRITE (6, *) ' at  z,ow' ' '=owW  (1-Rw2  )*,  z,  sigmawtp,  sigmawdp,  Rw 
d      IF(i.lt.iw) 

d     &WRITE(6, *)'au,,,,av,,, ,ow' ' *  * , sigmautp, sigmavtp, sigmawtp 
c  Get  random  normal  fluctuation  velocities,  utp  &  vtp,  at  time,  t-dt 

ruu  =  abs(pu/3.d0) 

rvv  =  abs(pv/3.d0) 

rww  =  abs(pw/3.d0) 
d     IF(i.lt.iw)  WRITE  (6, * ) ' ruu , rvv , rww ' , ruu, rvv, rww 

call  NORNG(  ruu,  pu  ) 

call  NORNG(  rvv,  pv  ) 

IF(  pw  .ge.  0.0  )wtp  =  upf actor* (pw* (sigmaw2  -  sigmawtp)  +  wbuoy) 

utp  =  sigmautp*pu 

vtp  =  sigmavtp*pv 

IF  (  noskew  .ne.  1  .and.  Linv  .It.  0.  )  GOTO  250 
call  NORNG(  rww,  pw  ) 
wtp  =  sigmawtp*pw  +  wbuoy 
d      IF(i.lt.iw)  WRITE  (6, * ) 'pu,pv,pw* ,pu,pv,pw 
d      IF(i.lt.iw)  WRITE  ( 6, *)' utp, vtp, wtp ', utp, vtp, wtp 
d      IF(i.lt.iw)  WRITE  (6f*)'w*"  =  ow1  "  *pw' , wtp, sigmawtp, pw 

GOTO  500 
2  50    CONTINUE 

d      IF(i.lt.iw)  WRITE  ( 6,  * )  ' pu,pv,pw' , pu, pv, pw 
d      IF(i.lt.iw)  WRITE  ( 6, *)' utp, vtp, wtp ', utp, vtp, wtp 

if  (  iMcNid  .ne.  1  )  GOTO  400 
c  modify  wtp  for  skewness  according  to  McNider  method 

iwp  =  0 
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iwm  =  0 
300    CONTINUE 

c  Get  random  normal  fluctuation  velocity,  wtp,  at  time,  t-dt 
call  NORNG(  rww,  pw  ) 
wtp  =  sigmawtp*pw 
IF  (  wtp  .ge.  0.  .and.  iwp  .eq. 
)  iwp  =  1 
.and.  iwm  .eq. 
)  iwm  =  1 


0. 
0. 
0. 
0 


or.  iwm  .eq.  0 


0  )  wp  =  wtp 
0  )  wm  =  wtp 
)  GOTO  300 


)S  =  0.1  -   0.2*zovrL**0.2 

)S  =  0.1  +  (0.6/(  0.68*(l.-15.*zovrL)**(-0.25) 


iwp,  iwm,  S  '  ,  iwp,  iwm,  S 


IF  (  wtp  .ge. 
IF  (  wtp  .le. 
IF  (  wtp  .le. 
IF  (  iwp  .eq. 
c  Get  S  function 

IF  (  zovrL  .gt.  0, 
IF  (  zovrL  .le.  0. 
&     -  1.8*zovrL  )  ) 
d      IF(i.lt.iw)  WRITE  (6,*) 
c  Get  alpha  and  eta 

alpha  =  -  0.028  -  0.6*abs(S) 
eta   =   0.54*abs(S) 
c  Get  wtp(t-dt)  !!!  Pielke  prints  this  on  p.  178  as  wdp(t-dt)  CHECK  THIS 
c  by  looking  at  McNider  paper!!!!! 

IF  (  zovrL  .le.  0.  )  wtp  =  alpha*wp  +  wm/alpha  -  eta 

IF  (  zovrl  .gt.  0.  )  wtp  =  wp/alpha  +  alpha*wm  +  eta 

d     IF(i.lt.iw)  WRITE  ( 6, *)' alpha, eta, wtp' , alpha, eta, wtp 

GOTO  500 
400    CONTINUE 

IF  (  iKamada  .ne. 


1  .and.  Linv  .It.  0.  )  GOTO  500 


c  Try  Kamada-Berkowicz  double  gaussian  skewness  approx  instead  of  McNider. 
c  This  approach  assumes  that  mean  absolute  vertical  fluctuation  velocity 
c  is   0.55wk,  vertically  averaged  over  the  BL  depth,  where  wk  generalizes 
c  w*  to  mixed  forced/free  convection.   We  let  the  up/downdraft  volume 
c  ratio  in  the  convective  BL  "  0.4/0.6,  with  mean  updraft  velocity 
c  "  1.08ow  and  the  mean  downdraft  velocity  "  0.87ow.   We  assume  that  the 
c  two  gaussian  velocity  distributions  are  displaced  from  zero  by  this 
c  amount  with  60%  of  the  distribution  in  the  downdraft  volume  and  40% 
c  in  the  updraft  volume.   We  assume  that  departures  from  Gaussian  are 
c  small  for  stable  cases.  (See  pps.  199-201,  Ch.  4  Lectures  on  Air 
c  Pollution  Modeling) 
c 

utp  =  pu*sigmaudp*sqrt (1  -  Ru**2) 
vtp  =  pv* sigma vdp*sqrt (1  -  Rv**2) 
wbuoy  =  Rw*wbuoy  -  (9 .801/Theta2 ) *Deltheta*dt 
sumwbuoy  =  sumwbuoy  +  wbuoy 
call  NORNG(  rww,  pw  ) 

IF  (  pw  .ge.  0.0  )  ikount  =  ikount  +  1 
d        IF(i.lt.iw)  WRITE  ( 6, *)' ikount , icount *, ikount , icount 
IF  (  ikount  .eq.  icount  )  THEN 
pw  =  -pw 
ikount  =  0 
ENDIF 
The  above  sets  an  (icount  +  1) / (2*icount )  chance  of  updraft  &  an 

(icount  -  1) / (2*icount )  chance  of  a  downdraft. 
E.g.,  if  icount  =  5,  then  updraft/downdraft  chances  are  60%/40%. 


c 
c 
c 
420 

c 


CONTINUE 
IF(  pw  .ge, 
IF(  pw  .le. 
IF(  pw  .ge. 


0.0  )wtp  =  upfactor*sigmawtp*pw  +  wbuoy 
0.0  )wtp  =  dnf actor*sigmawtp*pw  +  wbuoy 
0.0  ) 


i  wtp  =  upf actor*pw*sigmawdp*sqrt ( 1  -  Rw**2)  +  wbuoy 

IF(  pw  .le.  0.0  ) 

i  wtp  =  dnf actor*pw*sigmawdp*sqrt ( 1  -  Rw**2)  +  wbuoy 

IF(  pw  .le.  0.0  )  wtp  =  upf actor*pw*sigmawdp  +  wbuoy 
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d 

d 

d 

d 

500 

c 

c  Get 


IF(  pw  .le.  0.0  )  wtp  =  dnf actor*pw*sigmawdp  +  wbuoy 
sumw  =  sumw  +  wtp 

IF(i.lt.iw)WRITE(6,*)  'wtp,pw'  ,wtp,pw 
IF(iflag  .eq.  iw) 
&     WRITE ( 6, *) 'wbuoy, Rw, 50, sumw' ,  wbuoy, Rw,Deltheta, sumw 
IF(iflag  .eq.  ix)  WRITE ( 6, *)' sumwbuoy, sumw' ,  sumwbuoy, sumw 
CONTINUE 


fluctuation  velocities 

udp  =  udpold*Ru  +  utp 

vdp  =  vdpold*Rv  +  vtp 

wdp  =  wdpold*Rw  +  wtp 

IF  (ikamada  .eq.  1)  THEN 

udp  =  udpold  +  utp  -  (1.  -  Ru)*udpold 
vdp  =  vdpold  +  vtp  -  (1.  -  Rv)* vdpold 
wdp  =  wdpold  +  wtp  -  (1.  -  Rw)*wdpold 

ENDIF 
d      IF(i.lt.iw) 

d    &    WRITE  (6,*) 'at  z,wdp=wdpold*Rw+wtp' , z, wdp, wdpold, Rw, wtp 
c  Get  position 

x  =  x  +  udp*dt 

y  =  y  +  vdp*dt 

z  =  z  +  wdp*dt 


IF  (  z 

.It.  10*z0  )  THEN 

z  = 

abs(z)  +  10.*z0 

wdp 

=  abs(wdp) 

ELSEIF 

(  z  .gt.  zi  )  THEN 

z  = 

2.0*zi  -  z 

z  = 

MAX(z,  0.9*zi) 

wdp 

=  -abs(wdp) 

ENDIF 

udpold 

=  udp 

vdpold 

=  vdp 

wdpold 

=  wdp 

Tlw,  Rw,  Wz 


iflag  =  iflag  +  1 
IF  (iflag  .gt.  ix)  THEN 
d        WRITE  (6, 1110)time,Linv,Ri,udp, vdp,wdp,x,y,z 
d        WRITE  (6,1120)Km,  V,  ustar,  wpTpZ,  sigmaW2,  tke, 
iflag  =  1 
ENDIF 
*********************************************************************** 

c   debug  formats 
dlOOO  FORMAT  (lOx,  2gll.3) 
dlOlO  FORMAT  (lOx,  3fll.3) 
dl020  FORMAT  (lOx,  4fll.3) 
dl030  FORMAT  (lx,  i5  ) 
dl040  FORMAT  (lOx,  3gll.3,  i5) 
dl050  FORMAT  (lOx,  4fll.3,  i5) 
dl055  FORMAT  (lOx,  5fll.3,  i5) 
dl060   CONTINUE   . 
dlllO  FORMAT  (lx,'tm  1/L  Ri  u  v 
dll20  FORMAT  (lx,'Km  V  u*  Q  aw2 
RETURN 

END 
********************************************************************** 


w  x  y  z  ' ,f6.1,2f6.2,3f6.1,3f7.1) 
e  Tw  Rw  Wz  ' ,f6.1,5f6.2,f6.1,2f6.2) 


SUBROUTINE  MESOFLOW  (  z,  zi,  zO,  Linv,  windspd2 ,  ustar,  Km,  Ri, 
&      u,  v,  w,  Theta2,  dUtotdz,  dvdz,  wk,  D,  R,  wpTpZ,  BVF,  iw, 
&       BLS,  sigmaU2,  sigmaV2,  sigmaW2,  tke,  Wz,i,wk,  ustar, 
&      bdthetdz,  Thstar,  Deltheta) 
This  subroutine  computes  the  mesoscale  flow  parameters  used  as 
inputs  to  the  standard  mcnider  particle  calculations.   It  also  computes 
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c  some  alternate  values  for  the  au,  v,  w,  and  other  quantities. 

REAL  k,  LI,  InzovzO,  LovZi,  Km,  L,  lepsilon,  Ik 

DOUBLE  PRECISION  Linv 

DATA  g,  k,  zl,  f  /  9.801,  0.4,  2.0,  0.0001  / 
d      IF(i.lt.iw)WRITE(6,*) 'in  MESOFLOW  i  =  ',i 
c   Zero  out  some  variables 

u  =  0. 

v  =  0. 

w  =  0. 

dUdZ  =  0. 

dVdZ  =  0. 

wstar  =  0.0 

Wk  =  0. 

Ri  =  0. 

Km  =  0. 

tke  =  0. 

BVF  =  0. 

BLS  =  0. 
c   Initialize  some  values 

D  =  0.1 

R  =  0.2 

zOinv  =  l./zO 

InzovzO  =  alog( z*z0inv) 

IF  (  z  .It.  zO  )  z  =  z  +  zO 

zinv  =  l./z 

ziinv  =  l./zi 

L  =  l./Linv 

gzi  =  9.801*zi 

zlovL  =  zl*Linv 

Theta2  =  290.0 

zovL  =  z*Linv 
c  Get  psim  &  psih  for  surface  layer.   If  stable, 

psim  =  -5.*zlovL 

psih  =  psim 
d     IF (i.lt. iw)WRITE(6, *) 'psim=psih=-5zl/L' , psim, psih, zlovL 
c   If  unstable, 

IF  (  Linv  .It.  0.)  THEN 

phim  =  (1.  -  28.*zlovL)**(-0.25) 
d     IF(i.lt.iw)WRITE(6,*) 'phim=( l-28zl/L) ** (~k) ' ,phim,zlovL 
c        Get  psim  by  Kamada  (unpub) 

IF  (  zlovL  .It.  0  .and.  zlovL  .gt.  -0.01)  zlovL  =  -0.01 
d     IF(i.lt.iw)WRITE(6,*) ' just  before  psim  =  ' 

psim  =  (.0159651383  -  5 .4151107*zlovL) / 
&  (1.  -  3.59002231*zlovL  -0. 799168457*zlovL**2 ) 

d      IF(i.lt.iw)WRITE(6,*) ' psim= (a-czl/L) / ( l-bzl/L-d( zl/L) 2 ' ,psim,zlovL 
c        Also  from  Dyer  and  Bradley  (1982)  BLM  22,  3-19,  we  have  for  psi 
H 
d      IF(i.lt.iw)WRITE(6,*) ' just  before  psih  =  ' 

psih   =   2.*alog(    0.5*(1.    +   sqrtf    1.    -    14.*zlovL    ))     ) 
d  IF(i.lt.iw)WRITE(6,*)  'psih=2Ln(J{(  1+V  (l-14zl/L)  )  )  ■  ,psih,zlovL 

ENDIF 
c   Get  windspeed  drag  coefficient  and  friction  velocity 
d     IF(i.lt.iw)WRITE(6,*) ' just  before  sqrtcdn  =  ' 

sqrtCdn  =  0.4/(  alog( zl*z0inv)  -  psim) 
d      IF(i.lt.iw)WRITE(5,*) '/ Cdn=.4/ (Ln( zl/zO) -psim) ' , sqrtCdn, zl, zO, psim 

ustar  =  MAX(  sqrtCdn*windspd2,  0.01  ) 
d      IF(i.lt.iw)WRITE(6,*)  '  u*=MAX(\/  Cdn*V2  ,  0.  01)  '  ,  ustar,  sqrtCdn,  windspd2 
c 

c   For  dispersion,  get  vertical  profiles  of  relevant  variables,  using 
c   similarity  theory  from  ( Sorb j an  , Structure  of  the  Atmos  BL,  Ch.  4)  and 
c   (Panofsky  and  Dutton,  Atmos  Turbulence).   First  get  oft  used  variables 
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zovzi  =  z*ziinv 

zOovL  =  zO*Linv 

ziovL  =  zi*Linv 

Lovzi  =  L*ziinv 

ustar2  =  ustar*ustar 

ustar3  =  ustar*ustar2 

wpTpO  =  -0.25508*Linv*ustar3*Theta2 
d      IF(i.lt.iw)WRITE  ( 6, * ) ' wO=  =-. 255u**3©/L ' , wpTpO,ustar , Theta2 , L 

Tstar  =  -wpTpO/ustar 

IF  (  Linv  .It.  0.)  THEN 

Wstar  =  (gzi*wpTpO/Theta2)**. 33333 
d        IF(i.lt.iw)WRITE(6,*) 'w*=(gziw©/0) **l/3 ' ,wstar,gzi, wpTpO,Theta2 
c  Replace  w*  with  wk  for  mixed  forced/ free  convection  (Kamada, 

c  NPS-PH-92-007 .   Here,  add  only  surface  layer  shear  production. 

Ar  =  (1.  -  150. 0*z0ovL)**0. 333333  -  1 
d        IF (i. It. iw) WRITE  (6,*)'Ar=  (1  -  150z0/L) ** ( 1/3)  -  1' 
d         IF(i.lt.iw)WRITE  (6,*)Ar,  zOovL 
d      IF(i.lt.iw)WRITE(6,*) ' just  before  Wk  =  ' 

Wk  =  wstar*(l.  -  (  3.125  -  2.5*alog(Ar)  )*Lovzi  )**0. 33333 
d         IF(i.lt.iw)WRITE  (6, *) 'Wk=w* ( 1-3. 124-2 . 51n(Ar) )L/zi) ** ( 1/3 ) ' 
d         IF(i. It. iw)WRITE  ( 6, * ) wk, wstar , Ar ,Lovzi 

Thstar  =  -wpTpO/wk 
c   Decide  what  algos  to  use,  based  on  z/L. 
ENDIF 

IF  (  zovL  .It.  -0.5  )  icase  =  1 

IF  (  zovL  .It.  0.0  .and.  zovL  .ge.  -0.5  )  icase  =  2 

IF  (  zovL  .ge.  0.01  )  icase  =  3 

GOTO  (  330,  400,  430  )  icase 
330       CONTINUE 
********* 

c  Get  Km,  &  ou,v,w  for  forced/ free  convective  BL.   Ri  and  shear  not 

needed. 

********* 

c      Get  potential  temperature  at  height  z  above  -L/2 

R  =  abs(R) 

zovzil3  =  zovzi**. 33333 

zovzi23  =  zovzil3*zovzil3 

tpl  =  1.0  -  zovzi 

tp4  =  tpl  +  D 

tp413  =  tp4**. 33333 

tp423  =  tp4**. 66667 

Thstar  =  -wpTpO/wk 

tp23  =  tpl**. 666667 

R23  =  R**. 666667 

Oldtheta  =  Theta 

Theta  =  Theta2  +  Thstar* (  tp23/zovzil3  +  R23*zovzi23/tp413  ) 

IF  (  i.eg.  1)  Oldtheta  =  Theta 

Deltheta  =  Theta  -  Oldtheta 
d         IF(i.lt.iw)WRITE  ( 6, * ) ' 60, O, Old©' ,  Deltheta,  Theta,  Oldtheta 
c     temperature  and  buoyancy  fluxes  at  height,  z 

wpTpZ  =  wpTp0*(l.  -  1.2*zovzi) 

B  =  g/Theta2 

BwpTpZ  =  B*wpTpZ 
c      Get  vertical  velocity  variance  at  height  z  above  -L/2 

sigmaW2  =  wk**2* ( 1. l*zovzi23*tp23  +  R23*zovzi23*tp423 ) 
d      IF(i.lt.iw)WRITE  ( 6, * ) ' ow2=wk2 ( 1 . 1 ( z/zi) **2/3 ( 1-z/zi) **2/3  +  ' 
d      IF( i. It. iw) WRITE  (6,*)'R**(2/3)(z/zi)**(2/3) ( 1-z/zi+D) ** ( 1/3 ) ' 
d      IF(i. It. iw) WRITE  ( 6, * ) sigmaW2 , wk, zovzi, R,D 

sigmaW  =  sqrt ( sigmaW2 ) 
c      Get  TKE  dissipation  rate  at  height  z  above  -L/2,  using  similarity 

tmplO  =  Wk**3*ziinv 
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epsilonO=  0.75*tmplO 

phiEp  =  0.75*tpl  +  R*zovzi 

epsilon  =  phiEp*tmplO 
d      IF(i. It. iw) WRITE  (6,*) 'e  =  *e*wk**3/zi ', epsilon, phiEp,wk, zi 
c      Get  ou2,av2 

wpsqovE  =  0.333*(  (2.*epsilon0  -  epsilon) /epsilonO  + 
&  BwpTpZ /epsilon  ) 

d      IF(i.lt.iw)WRITE(6,*) 'w'2/e=(2e0-e)/3e0  +^'0'  z/e  '  ,wpsqovE 
d      IF(i.lt.iw)WRITE(6,*) ■ eO,  fiw'O'z*,  epsilonO,  BwpTpZ 

sigmaU2  =  0. 5555*sigmaW2* (  2 . /wpsqovE  -  1.) 
c         sigmaU2  =  MAX(  sigmaU2,  .33333  ) 

sigmaV2  =  0.8*sigmaU2 
d      IF(i. It. iw) WRITE  ( 6, * ) sigmaU2 , wk, zovzi, sigmW2 
c      The  following  tke  parameterization  is  OK  for  low  shear 

tke  =  0.5*(  sigmaU2  +  sigmaV2  +  sigmaW2  ) 
d     IF(i.lt.iw)WRITE  (6,*)*e  =  ow2/2  +  3ou2/2 ' , tke, sigmaW2 , sigmaU2 

tkel2inv  =  l./sqrt(tke) 

tke32inv  =  tkel2inv**3 
c     Get  dissipation  &  mixing  length  scales  &  vertical  diffusion 
c      coefficient 

lepsilon  =  l./(  zinv  +  1. 4*epsilon*tke32inv  ) 
d     IF(i.lt.iw)WRITE(6,*)  '  le  =  l/(l/z  +1.4e</e  )', lepsilon, z, epsilon, tke 

Ik  =  MIN  (  z,  3.0*lepsilon*sigmaW2/tke  ) 
d     IF(i.lt.iw)WRITE  ( 6, * ) ' Ik  =  MIN( 31e*ow2/tke,  z) ', Ik, lepsilon, tke, z 
c     The  mixing  length  scale  and  diffusion  coefficients  are  Kamada 
c     variants 

Km  =  0.44*lk/tkel2inv 
d      IF(i. It. iw) WRITE  (6,*) 'Km  =  0. 44*lk*sqrt ( tke) ' ,Km, Ik, tke 
c  Get  temperature  and  buoyancy  gradients 

c     dthetadz  =  . 6*Thstar*ziinv* (tp23/zovzi23**2  -2. *R23*zovzi23/tp423 ) 
dthetadz  =  . 6*Thstar*ziinv* (tp23/zovzil3    -2 . *R23*zovzi23/tp423) 
Bdthetdz  =  B*dthetadz 
GOTO  460 

400    CONTINUE 
********** 

c   get  Km  for  unstable  surface  layer.    Shear  and  Ri  not  needed  for 

unstable. 

********** 

Km  =  k*ustar*z/phim 

tkeusl  =  (5.6*Km*zinv)**2 

Ri  =  zovL 

phil  =  (  12.  -  0.5*ziovL  )**. 333333 

phi2  =  phil 

phi3  =  (  1.  -  14.*zovL  )**.25 

phiEp=  (1.  +  .75*(-zovL)**.6667)**1.5 

epsilonO  =  ustar3*phiEp*zinv 

GOTO  440 

430    CONTINUE 
********** 

c   get  Km,  Ri,  and  other  parameters  for  neutral  to  stable  BL 
********** 

h  =  zi 

IF  (  zi  .It.  -98.  )  h  =  sqrt(k*ustar*Linv/f ) 
d        IF(i.lt.iw)WRITE(6,*) 'h=zi,or  / (ku*/(fL) ) ' ,h,zi, f 

hinv  =  l./h 

zovH  =  z*hinv 

zsqovHL  =  zovL*zovH 

zsqovH2  =  zovH*zovH 

alpha  =  2.0  -  10.0*Linv 

beta   =  3.0  -  20.0*Linv 
d         IF(i.lt.iw)WRITE(6,*)  'a,/3,z/H'  , alpha,  beta,  zovH 
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tp8  =  1.  -  zovH 

tp8alpha  =  tp8**alpha 
d         IF(i.lt.iw)WRITE(6,*) ' (1-z/H) , (l-z/H)**a' ,  tp8,tp8alpha 

tp8al2   =  sqrt (tp8alpha) 

tp8beta  =  tp8**beta 
d         IF(i.lt.iw)WRITE(6,*)  '  ( 1-z/H)  **a/2,  ( 1-z/H)  **J3'  ,  tp8al2  ,  tp8beta 

phil  =  2.2*tp8al2 

phi2  =  phil 

phi3  =  1.6*tp8al2 
c      local  friction  velocity 

ul  =  ustar*sqrt (tp8alpha) 
d        IF(i.lt.iw)WRITE(6,*) *ul=u* ( 1-z/H) ** (a/2 ) ' ,ul,ustar, zovH 
c      local  Obukhov  length 

LI  =  L*tp8**(1.5*alpha  -  beta) 
d         IF(i.lt.iw)WRITE(6,*)  'L1=L(  1-z/H)  **  (3a/2-J3)  '  ,L1,L,  zovH,  alpha, 
d    &  beta 

c      total  shear 
c        dUtotdz  =  4.7*ul/(k*Ll) 

dUtotdz  =  2.5*ustar*zinv*(l.+4.7*zovL)**(0.5*alpha)*tp8al2 
d         IF(i.lt.iw)WRITE(6,*) 'dU/dz=4. 7ul/ (kLl) ' , dUtotdz, ul, k,Ll 
c     temperature  and  buoyancy  fluxes  at  height,  z 

wpTpZ  =  wpTpO*(l.  -  zovzi)**beta 
c      Brunt  Vaisala  frequency  at  z 

tp9  =  (1.  +  3.7*zovL)*zinv 
d        IF(i.lt.iw)WRITE(6,*) • tp9=( l+3.7z/L) /z ' , tp9, zovL, zinv 

BVF  =  4.3*ustar*Tstar**2*tp9*tp8beta*tp8alpha 
d         IF(i.lt.iw)WRITE(6,*)  '  BVF=4.  3u*0*2  ( 1+3  .  7z/L)  ( 1-z/H)  **  (a+13) /z  *  , 
d    &   BVF,ustar,Tstar, zovH, zovL, z 
c      TKE  dissipation  rate  at  z 

phiEp  =  3.6*tp9*tp8beta*zinv 

epsilon  =  phiEp*ustar3 
c      Richardson  number  at  z 

tplO  =  l./(l.  +  5.*zovL  ) 

Ri  =  zovL*tplO 
c      momentum/heat  dif fusivities  at  z 

Km  =  k*ustar*z*tp8*tpl0 
d         IF(i.lt.iw)WRITE(6,*) ' Km=ku*z ( 1-z/H) /( l+5z/L) ' ,Km, k,ustar, z, 
d     &  zovH,zovL 

***************************** 

440       CONTINUE 

c  Get  ou,v,w,  tke,  total  horizontal  velocity,  &  temperature 

sigmaU  =  ustar*phil 

sigmaU2=  sigmaU*sigmaU 

sigmaV  =  0.894*sigmaU 

sigmaV2=  0.5*sigmaU2 

sigmaW  =  ustar*phi3 

sigmaW2  =  sigmaW**2 

tke  =  0.5*(sigmaU2  +  sigmaV2  +  sigmaW2) 
d         IF(i.lt.iw)WRITE(6,*) ' tke=oui2/2 ' , tke, sigmaU2, sigmaV2 , sigmaW2 
c        IF  (  icase  .eq.  2  )  tke  =  tkeusl 

Oldtheta  =  Theta 

Theta  =   Theta2  +  2.5*Tstar*(  InzovzO  -  psih  ) 

IF  (  i.  eq.  1)  Oldtheta  =  Theta 

Deltheta  =  Theta  -  Oldtheta 
d        IF(i.lt.iw)WRITE  ( 6, *) ' 5O,O,Old0' ,  Deltheta,  Theta,  Oldtheta 
460       CONTINUE 

Vtotal  =  2. 5*ustar* (InzovzO  -  psim) 
d         IF(i.lt.iw)WRITE(6,*) ' V=2 . 5u* (Ln( z/zO) -psim) ) ' , Vtotal, ustar, 
d     &  InzovzO, psim 

u  =  Vtotal 
d         IF(i.lt.iw)WRITE(6,*) '0=02+2 . 50* (Ln( z/z0)-psih) ' , theta, theta2 , 
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d    &  tstar, lnzovzO,psih 

********************************************************************** 

c   debug  formats 

dlOOO  FORMAT  (lOx,  2gll.3) 

dlOlO  FORMAT  (lOx,  3fll.3) 

dl020  FORMAT  (lOx,  4fll.3) 

dl030  FORMAT  (lx,  i5  ) 

dl040  FORMAT  (lOx,  3gll.3,  i5) 

dl050  FORMAT  (lOx,  4fll.3,  i5) 

dl055  FORMAT  (lOx,  5fll.3,  i5) 

dl060   CONTINUE 

RETURN 

END 
********************************************************************** 

SUBROUTINE  NORNG(r,p) 
c 

c   This  subroutine  generates  a  sequence  of  numbers  normally  and  randomly 
c  distributed  over  the  interval  -3  to  3  from  uniformly  distributed  random 
c   numbers,  by  the  method  of  linear  approximation  to  the  inverse  of  the 
c   accumulative  normal  distribution  function 
c 

DOUBLE  PRECISION  r,p,  y(6),  x(6),  s(5) 

DATA  y/O.dO,  0.0228d0,  0.0668d0,  0.1357d0,  0.2743d0,  0.5d0  /, 
&      x/  -3.01d0,  -2.0d0,  -1.5d0,  -l.OdO,  -0.6d0,  O.OdO  /, 
&  3/    43.8596d0,  11.3636d0,  7.25689d0,  2.891352d0,  2.65887d0  / 

CALL  STRNUM(r) 

p  =  r 

i  =  1 

IF  (  p  .gt.  0.5d0  )  p  =  l.OdO  -  r2      IF  (  p  .It.  y(i+l)  )  GOTO  8 

i  =  i  +  1 

GOTO  2 
8     P  =  (  (  P  -  y(i)  )*s(i)  +  x(i)  ) 

IF  (  r  .ge.  0.5d0  )  p  =  -p 
cd    WRITE  (6,*)' random  numbers,  r  &  p',r,p 

RETURN 

END 
************************************************************************ 

SUBROUTINE  STRNUM(r) 
c 

c  This  subroutine  generates  a  sequence  of  numbers  which  are  randomly  and 
c   uniformly  distributed  over  the  unit  interval, 
c 

DOUBLE  PRECISION  pi,  r 

bb  =  l.dO 

pi  =  r*317.d0 
c      VfRITE  (6,*)  'pi'  , pi 

r   =  abs(  pi  -  (  INT(pl/bb)*bb)  ) 

RETURN 

END 
************************************************************************ 

FUNCTION  MAX  (a,  b) 
max  =  b 

IF  (  a  .gt.  b)  max  =  a 
END 
************************************************************************ 

FUNCTION  MIN  (a,  b) 
min  =  b 

IF  (  a  .It.  b)  max  =  a 
END 
************************************************************************ 
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